简介这份资源面向电力系统状态估计与网络安全方向的研究生、科研人员及工程技术人员聚焦虚假数据注入攻击的防御问题提供基于鲁棒广义极大似然GM估计器的完整MATLAB实现方案。资源包共13个文件以10个m脚本为核心配合1份pdf说明、1份docx文档及1个txt许可文件整体约159KB涵盖GM估计器主流程、Givens旋转数值稳定化、零注入处理与变压器抽头联合估计等关键模块。已有1104人学习下载说明其在电力监控与网络攻防领域具有较高参考价值。读者可据此复现投影统计鲁棒估计的完整流程理解坏数据、坏杠杆点与恶意注入攻击下的防御机理并借助脚本中的修正因子、稀疏矩阵与测试对比代码快速搭建仿真环境、验证算法在高斯及厚尾噪声下的统计效率为课题研究与工程应用提供可直接运行的基础代码与排错思路。1. 从一次变电站遥测跳变说起鲁棒电力系统状态估计器到底在防什么凌晨两点某地区调度中心的 SCADA 画面上一条 220kV 线路的有功潮流在 40 秒内从 180MW 跳到 420MW又跳回来量测通道自检全绿通信误码率为零。值班员第一反应是 CT 饱和或者通道抖动但核对相邻变电站的对应量测后发现两侧功率不守恒——这不是设备问题是有人往量测里塞了假数据。这类场景对应的技术名词就是虚假数据注入攻击False Data Injection AttackFDIA而用来扛住它的核心组件就是鲁棒电力系统状态估计器。传统加权最小二乘WLS状态估计有个致命前提量测噪声服从零均值高斯分布且坏数据可以被残差检测出来。FDIA 的高明之处在于攻击者如果掌握网络拓扑和支路参数可以构造一个注入向量让状态估计的残差几乎不变但状态变量被系统性偏移。换句话说坏数据检测这个黑匣子被绕过去了。鲁棒状态估计器的思路不是去猜攻击者怎么构造向量而是换一套估计准则让少量被篡改的量测无法主导整个解。这套方法适合三类人做电网调度自动化的工程师、研究信息物理系统安全的团队、以及需要给状态估计模块做加固的二次开发人员。下面从选型、实现到踩坑把这条路走一遍。2. 鲁棒状态估计器的选型为什么不是简单换个损失函数2.1 WLS、WLAV、GM 估计器的本质差别很多人第一次接触鲁棒估计会以为把 WLS 的目标函数从平方换成绝对值就行。方向对但不够。WLS 的目标是最小化加权残差平方和它对大残差极度敏感——一个被篡改的量测残差翻三倍代价函数贡献翻九倍解会被它拽着走。加权最小绝对值WLAV把平方换成绝对值对大残差的惩罚从二次降为一次抗差能力立刻上一个台阶。但 WLAV 在残差接近零时不可导数值求解要用线性规划收敛速度慢而且当坏数据比例超过某个阈值时会突然失稳。广义极大似然GM估计器走的是另一条路用一个有界的影响函数influence function来压制大残差。常见的有 Huber、Tukey 双权、Hampel 三段式。Huber 在残差小于阈值时保持二次超过阈值后转为线性兼顾了正常量测的效率和异常量测的鲁棒性。Tukey 更激进超过阈值直接给零权重相当于把可疑量测踢出估计。选哪个取决于你能容忍多少量测被误杀。估计器目标函数抗差机制求解方式适用场景WLSΣ wᵢrᵢ²无正规方程/牛顿法无攻击、噪声干净WLAVΣ wᵢ|rᵢ|线性惩罚线性规划坏数据比例10%Huber GM分段二次/线性阈值截断迭代重加权最小二乘坏数据比例10%~30%Tukey GM有界红降函数零权重剔除迭代重加权坏数据比例高但稀疏我一般会先上 Huber因为它的阈值参数有明确的统计解释阈值取 1.5 倍量测标准差时正常量测被误判的概率约 13%但权重只降到 0.7 左右不会直接丢掉。如果攻击者注入的量测比例超过 20%再考虑 Tukey 或者引入投影统计量做初值筛选。2.2 投影统计量给状态估计器装一个预筛层GM 估计器有个隐患迭代重加权依赖初值如果初值被坏数据带偏后面再鲁棒也拉不回来。投影统计量Projection Statistics就是解决初值问题的。它的思路是把每个量测向量往多个方向上投影正常量测的投影应该聚集在某个范围偏离太远的直接标记为可疑。这一步不求解状态变量只做量测空间的离群检测计算量小可以放在状态估计之前。具体做法是对量测矩阵的每一列对应一个量测计算它在所有可能投影方向上的中位数和绝对偏差得到一个稳健的马氏距离。超过卡方分布阈值的量测进入可疑集在后续 GM 迭代中给它们更低的初始权重。这一步相当于给状态估计器加了一个预检门把明显离谱的量测挡在外面避免它们污染初值。注意投影统计量的计算复杂度是 O(m²)m 是量测数。对于几千个量测的区域电网这一步可能比状态估计本身还慢。常见做法是只对残差最大的前 20% 量测做投影统计或者用随机投影降维。2.3 量测冗余度鲁棒估计器的生命线再鲁棒的估计器也怕量测不够。如果某个节点的注入功率只有一个量测攻击者改它估计器没有任何交叉验证的依据。电力系统状态估计的可观测性分析里有个关键指标叫冗余度等于量测数除以状态变量数。冗余度低于 1.5 时鲁棒估计器的效果会急剧下降因为坏数据检测的自由度不够。我见过一个 14 节点系统量测配置只覆盖了 80% 的支路功率冗余度 1.2。在这种配置下Huber 估计器和 WLS 的差别不到 5%攻击者只要改两个关键量测就能把状态拉偏。后来补了 PMU 的量测冗余度提到 2.1同样的攻击场景下 Huber 估计器的状态偏差从 12% 降到 3% 以内。所以做鲁棒估计之前先算冗余度低于 1.8 的话优先补量测而不是调算法。3. 用 Python 跑通一个最小鲁棒状态估计器3.1 构造 IEEE 14 节点算例与量测向量先搭一个能复现的算例。用 pandapower 建 IEEE 14 节点模型生成潮流真值再按真值加高斯噪声造量测最后注入虚假数据。这一步的关键是攻击向量要满足 FDIA 的构造条件即攻击后的量测残差与攻击前几乎一致否则随便加个噪声都能被检测出来测不出鲁棒估计器的真实能力。import numpy as np import pandapower as pp import pandapower.networks as pn # 建 IEEE 14 节点模型并跑潮流 net pn.case14() pp.runpp(net) # 提取真值节点电压幅值、相角支路功率 V_true net.res_bus.vm_pu.values theta_true np.deg2rad(net.res_bus.va_degree.values) P_branch_true net.res_line.p_from_mw.values Q_branch_true net.res_line.q_from_mw.values # 组装量测向量 z h(x) e这里简化为直接用量测函数 # 实际工程中 h(x) 是非线性潮流方程这里用真值加噪声模拟 np.random.seed(42) sigma 0.01 # 量测噪声标准差 z_voltage V_true np.random.normal(0, sigma, len(V_true)) z_power P_branch_true np.random.normal(0, sigma * 100, len(P_branch_true)) # 构造 FDIA 攻击向量攻击者篡改 3 号和 8 号节点的注入功率量测 # 攻击量 a 满足 a H * cc 是状态偏移向量H 是量测雅可比矩阵 # 这里简化处理直接在量测上加一个与拓扑相关的偏移 attack_idx [2, 7] # 对应节点 3 和 8 z_power_attacked z_power.copy() z_power_attacked[attack_idx] np.array([15.0, -12.0]) # 注入虚假功率偏移 print(f攻击前量测均值: {z_power.mean():.2f}) print(f攻击后量测均值: {z_power_attacked.mean():.2f})这段代码做了三件事跑潮流拿真值、加噪声造量测、在指定量测上注入偏移。参数sigma控制噪声水平实际工程中功率量测的噪声标准差通常在 1%~2% 额定值这里用sigma * 100是因为功率基准是 100MW。攻击偏移量 15MW 和 -12MW 是随手设的真实攻击者会按a Hc构造让残差不变但这里为了演示鲁棒估计器的压制效果直接用固定偏移就够了。3.2 Huber 估计器的迭代重加权实现Huber 估计器的核心是迭代重加权最小二乘IRLS。每一轮用当前残差算权重残差大的量测权重低然后解一次 WLS更新状态再算残差直到收敛。下面是一个简化版实现状态变量只取电压幅值和相角量测函数用线性化近似。def huber_weight(residual, delta1.5): Huber 权重函数残差小于 delta 时权重为 1超过时按 delta/|r| 衰减 abs_r np.abs(residual) weights np.ones_like(abs_r) mask abs_r delta weights[mask] delta / abs_r[mask] return weights def robust_state_estimation(z, H, x0, max_iter20, tol1e-6, delta1.5): z: 量测向量 (m,) H: 量测雅可比矩阵 (m, n) x0: 状态初值 (n,) delta: Huber 阈值通常取 1.5 倍量测标准差 x x0.copy() for it in range(max_iter): # 计算残差 r z - H x # 算 Huber 权重 w huber_weight(r, delta) # 加权最小二乘解x (H^T W H)^(-1) H^T W z W np.diag(w) HtWH H.T W H HtWz H.T W z x_new np.linalg.solve(HtWH, HtWz) # 收敛判断 if np.linalg.norm(x_new - x) tol: print(f收敛于第 {it1} 次迭代) break x x_new return x, w # 构造简化的量测雅可比矩阵实际应用需按潮流方程求偏导 # 这里用随机矩阵模拟仅演示算法流程 m, n 20, 10 H np.random.randn(m, n) x_true np.random.randn(n) z_clean H x_true np.random.normal(0, 0.01, m) z_attack z_clean.copy() z_attack[2] 0.5 # 注入攻击 z_attack[7] - 0.4 x0 np.zeros(n) x_est_clean, w_clean robust_state_estimation(z_clean, H, x0) x_est_attack, w_attack robust_state_estimation(z_attack, H, x0) print(f干净数据状态误差: {np.linalg.norm(x_est_clean - x_true):.4f}) print(f攻击数据状态误差: {np.linalg.norm(x_est_attack - x_true):.4f}) print(f被攻击量测的权重: {w_attack[2]:.3f}, {w_attack[7]:.3f})这段代码里huber_weight是权重函数delta是阈值取 1.5 倍量测标准差是经验值。robust_state_estimation做 IRLS 迭代每次用当前残差更新权重再解加权最小二乘。关键参数max_iter控制最大迭代次数tol是收敛容差。运行后你会看到干净数据下状态误差很小攻击数据下误差被压制而且被攻击量测的权重明显低于 1。这就是鲁棒估计器在起作用——它没有去识别哪个量测被攻击而是通过降权让攻击量测无法主导解。3.3 用残差协方差做攻击检测的辅助判据鲁棒估计器本身不输出有没有攻击的结论它只是让估计结果更稳。如果你需要报警还得加一个检测环节。常用的是归一化残差检验算每个量测的残差除以其标准差超过阈值就报警。但 FDIA 的残差可能很小所以更可靠的是用鲁棒估计器的权重分布——如果大量量测权重同时下降说明系统里存在系统性偏差而不是单个坏数据。def attack_detection(w, threshold0.5, ratio0.3): w: 鲁棒估计器输出的权重向量 threshold: 权重低于此值视为可疑 ratio: 可疑量测比例超过此值触发报警 suspicious np.sum(w threshold) suspicious_ratio suspicious / len(w) if suspicious_ratio ratio: return True, suspicious_ratio return False, suspicious_ratio # 用上面的权重做检测 is_attack_clean, ratio_clean attack_detection(w_clean) is_attack, ratio_attack attack_detection(w_attack) print(f干净数据可疑比例: {ratio_clean:.2%}, 报警: {is_attack_clean}) print(f攻击数据可疑比例: {ratio_attack:.2%}, 报警: {is_attack})这个检测逻辑很简单统计权重低于 0.5 的量测比例超过 30% 就报警。参数threshold和ratio需要根据实际系统的量测冗余度和噪声水平调。冗余度高的系统可以放宽ratio因为正常量测多少数被降权不影响比例。冗余度低的系统要收紧否则容易漏报。4. 避坑与排查鲁棒状态估计器落地时的五个血泪教训4.1 现象估计结果震荡不收敛迭代 50 次还在跳原因Huber 阈值delta设得太小正常量测也被降权权重矩阵每轮剧烈变化IRLS 在解附近来回震荡。或者量测雅可比矩阵H的条件数太大加权后更病态。解决先把delta调到 2.0~2.5 倍量测标准差观察收敛曲线。如果还震荡检查H矩阵的条件数超过 1e6 的话需要做量测筛选或加正则化项。我一般会在HtWH上加一个小的对角项1e-6 * I相当于岭回归能显著改善数值稳定性。4.2 现象攻击量测的权重没降下来估计结果还是被带偏原因攻击者构造的虚假数据与正常量测的残差分布很接近Huber 权重函数在阈值附近区分度不够。或者攻击量测的数量超过了鲁棒估计器的崩溃点breakdown pointHuber 的崩溃点约 50%但实际有效范围通常只有 30%。解决换 Tukey 双权函数它的红降特性对接近阈值的残差更敏感。或者引入投影统计量做预筛把可疑量测在迭代前就标记出来给它们更低的初始权重。如果攻击量测比例确实超过 30%单靠鲁棒估计器不够需要结合 PMU 的动态量测做交叉验证。4.3 现象投影统计量计算太慢实时性达不到要求原因投影统计量要对每个量测计算所有投影方向的中位数和 MAD复杂度 O(m²)m 是量测数。区域电网 m 可能上千单次计算就超过状态估计本身的时间。解决只对残差最大的前 20% 量测做投影统计其余量测直接给正常权重。或者用随机投影代替全方向投影随机选 50~100 个方向精度损失很小但速度提升一个数量级。另一个做法是把投影统计量放在状态估计之前做一次后续迭代不再重复计算。4.4 现象量测冗余度不足时鲁棒估计器和 WLS 结果几乎一样原因冗余度低于 1.5 时坏数据检测的自由度不够鲁棒估计器的权重调整空间被压缩。攻击者只要改少数关键量测就能同时骗过 WLS 和鲁棒估计器。解决优先补量测尤其是 PMU 的电压相角量测它对状态估计的可观测性贡献最大。如果补不了量测退而求其次用历史数据做时序一致性检验——攻击者可以改单点量测但很难同时改多个时间断面的量测而保持时序连贯。把时序残差也纳入权重计算能部分弥补冗余度不足。4.5 现象攻击检测误报率高正常操作也被报警原因检测阈值ratio设得太低或者系统本身存在量测偏差比如 CT 慢漂移导致正常量测的权重也偏低。另外如果系统里有大量零注入节点这些节点的量测权重天然不稳定容易触发误报。解决先做一轮无攻击场景的基线测试统计正常情况下的可疑量测比例把ratio设成基线的 2~3 倍。对零注入节点单独处理不纳入可疑比例统计。如果 CT 漂移是已知问题在状态估计之前先做量测校准别让鲁棒估计器去扛这个锅。5. 从离线验证到在线部署一个可复用的验证套路鲁棒状态估计器写完只是第一步怎么证明它在真实攻击下有效才是决定要不要投入的关键。我一般会走三步验证离线注入测试、半实物仿真、现场试运行。离线测试用历史量测数据人为注入不同比例的 FDIA看状态偏差和检测率。半实物仿真用 RTDS 或者 RT-LAB 接真实 PMU验证通信延迟和量测丢包对鲁棒估计器的影响。现场试运行先旁路运行不接入闭环控制只记录估计结果和报警日志跑两周再评估。下面是一个离线验证的脚本框架用蒙特卡洛跑 100 次不同攻击场景统计状态误差和检测率。def monte_carlo_validation(n_trials100, attack_ratio0.2): n_trials: 蒙特卡洛次数 attack_ratio: 被攻击量测的比例 errors_wls [] errors_huber [] detection_rates [] for trial in range(n_trials): # 每次重新生成量测和攻击 m, n 30, 12 H np.random.randn(m, n) x_true np.random.randn(n) z H x_true np.random.normal(0, 0.01, m) # 随机选 attack_ratio 比例的量测注入攻击 n_attack int(m * attack_ratio) attack_idx np.random.choice(m, n_attack, replaceFalse) z_attack z.copy() z_attack[attack_idx] np.random.normal(0, 0.5, n_attack) # WLS 估计 x_wls np.linalg.lstsq(H, z_attack, rcondNone)[0] errors_wls.append(np.linalg.norm(x_wls - x_true)) # Huber 估计 x_huber, w robust_state_estimation(z_attack, H, np.zeros(n)) errors_huber.append(np.linalg.norm(x_huber - x_true)) # 检测 is_attack, _ attack_detection(w) detection_rates.append(1 if is_attack else 0) print(fWLS 平均状态误差: {np.mean(errors_wls):.4f}) print(fHuber 平均状态误差: {np.mean(errors_huber):.4f}) print(f攻击检测率: {np.mean(detection_rates):.2%}) print(f误差降低幅度: {(1 - np.mean(errors_huber)/np.mean(errors_wls)):.2%}) monte_carlo_validation(n_trials100, attack_ratio0.2)这个脚本跑 100 次每次随机选 20% 的量测注入攻击对比 WLS 和 Huber 的状态误差。参数attack_ratio可以调从 0.1 到 0.4 各跑一遍看鲁棒估计器的误差降低幅度怎么变化。如果attack_ratio超过 0.3 后误差降低幅度骤降说明这个配置下的崩溃点到了需要补量测或者换更强的鲁棒估计器。验证通过后在线部署还有几个工程细节要注意。第一状态估计的周期通常是 5~15 秒鲁棒估计器的迭代次数要控制在这个时间窗内max_iter别超过 10。第二权重矩阵W的存储和计算要优化用稀疏矩阵别用稠密np.diag。第三报警日志要记录每次迭代的权重分布方便事后回溯——攻击者可能慢慢调大量测偏移单次看不出来但权重分布的趋势会暴露问题。我自己的习惯是每次现场试运行前先用历史数据跑一遍离线验证把delta、ratio、max_iter这三个参数记在配置文件里别硬编码在代码里。现场环境一变量测噪声水平可能差一倍参数不调的话鲁棒估计器要么不收敛要么误报率飙升。这套东西没有一劳永逸的参数只有不断根据现场数据微调的习惯。希望帮到你。本文还有配套的精品资源点击获取
简介这份资源面向电力系统状态估计与网络安全方向的研究生、科研人员及工程技术人员聚焦虚假数据注入攻击的防御问题提供基于鲁棒广义极大似然GM估计器的完整MATLAB实现方案。资源包共13个文件以10个m脚本为核心配合1份pdf说明、1份docx文档及1个txt许可文件整体约159KB涵盖GM估计器主流程、Givens旋转数值稳定化、零注入处理与变压器抽头联合估计等关键模块。已有1104人学习下载说明其在电力监控与网络攻防领域具有较高参考价值。读者可据此复现投影统计鲁棒估计的完整流程理解坏数据、坏杠杆点与恶意注入攻击下的防御机理并借助脚本中的修正因子、稀疏矩阵与测试对比代码快速搭建仿真环境、验证算法在高斯及厚尾噪声下的统计效率为课题研究与工程应用提供可直接运行的基础代码与排错思路。1. 从一次变电站遥测跳变说起鲁棒电力系统状态估计器到底在防什么凌晨两点某地区调度中心的 SCADA 画面上一条 220kV 线路的有功潮流在 40 秒内从 180MW 跳到 420MW又跳回来量测通道自检全绿通信误码率为零。值班员第一反应是 CT 饱和或者通道抖动但核对相邻变电站的对应量测后发现两侧功率不守恒——这不是设备问题是有人往量测里塞了假数据。这类场景对应的技术名词就是虚假数据注入攻击False Data Injection AttackFDIA而用来扛住它的核心组件就是鲁棒电力系统状态估计器。传统加权最小二乘WLS状态估计有个致命前提量测噪声服从零均值高斯分布且坏数据可以被残差检测出来。FDIA 的高明之处在于攻击者如果掌握网络拓扑和支路参数可以构造一个注入向量让状态估计的残差几乎不变但状态变量被系统性偏移。换句话说坏数据检测这个黑匣子被绕过去了。鲁棒状态估计器的思路不是去猜攻击者怎么构造向量而是换一套估计准则让少量被篡改的量测无法主导整个解。这套方法适合三类人做电网调度自动化的工程师、研究信息物理系统安全的团队、以及需要给状态估计模块做加固的二次开发人员。下面从选型、实现到踩坑把这条路走一遍。2. 鲁棒状态估计器的选型为什么不是简单换个损失函数2.1 WLS、WLAV、GM 估计器的本质差别很多人第一次接触鲁棒估计会以为把 WLS 的目标函数从平方换成绝对值就行。方向对但不够。WLS 的目标是最小化加权残差平方和它对大残差极度敏感——一个被篡改的量测残差翻三倍代价函数贡献翻九倍解会被它拽着走。加权最小绝对值WLAV把平方换成绝对值对大残差的惩罚从二次降为一次抗差能力立刻上一个台阶。但 WLAV 在残差接近零时不可导数值求解要用线性规划收敛速度慢而且当坏数据比例超过某个阈值时会突然失稳。广义极大似然GM估计器走的是另一条路用一个有界的影响函数influence function来压制大残差。常见的有 Huber、Tukey 双权、Hampel 三段式。Huber 在残差小于阈值时保持二次超过阈值后转为线性兼顾了正常量测的效率和异常量测的鲁棒性。Tukey 更激进超过阈值直接给零权重相当于把可疑量测踢出估计。选哪个取决于你能容忍多少量测被误杀。估计器目标函数抗差机制求解方式适用场景WLSΣ wᵢrᵢ²无正规方程/牛顿法无攻击、噪声干净WLAVΣ wᵢ|rᵢ|线性惩罚线性规划坏数据比例10%Huber GM分段二次/线性阈值截断迭代重加权最小二乘坏数据比例10%~30%Tukey GM有界红降函数零权重剔除迭代重加权坏数据比例高但稀疏我一般会先上 Huber因为它的阈值参数有明确的统计解释阈值取 1.5 倍量测标准差时正常量测被误判的概率约 13%但权重只降到 0.7 左右不会直接丢掉。如果攻击者注入的量测比例超过 20%再考虑 Tukey 或者引入投影统计量做初值筛选。2.2 投影统计量给状态估计器装一个预筛层GM 估计器有个隐患迭代重加权依赖初值如果初值被坏数据带偏后面再鲁棒也拉不回来。投影统计量Projection Statistics就是解决初值问题的。它的思路是把每个量测向量往多个方向上投影正常量测的投影应该聚集在某个范围偏离太远的直接标记为可疑。这一步不求解状态变量只做量测空间的离群检测计算量小可以放在状态估计之前。具体做法是对量测矩阵的每一列对应一个量测计算它在所有可能投影方向上的中位数和绝对偏差得到一个稳健的马氏距离。超过卡方分布阈值的量测进入可疑集在后续 GM 迭代中给它们更低的初始权重。这一步相当于给状态估计器加了一个预检门把明显离谱的量测挡在外面避免它们污染初值。注意投影统计量的计算复杂度是 O(m²)m 是量测数。对于几千个量测的区域电网这一步可能比状态估计本身还慢。常见做法是只对残差最大的前 20% 量测做投影统计或者用随机投影降维。2.3 量测冗余度鲁棒估计器的生命线再鲁棒的估计器也怕量测不够。如果某个节点的注入功率只有一个量测攻击者改它估计器没有任何交叉验证的依据。电力系统状态估计的可观测性分析里有个关键指标叫冗余度等于量测数除以状态变量数。冗余度低于 1.5 时鲁棒估计器的效果会急剧下降因为坏数据检测的自由度不够。我见过一个 14 节点系统量测配置只覆盖了 80% 的支路功率冗余度 1.2。在这种配置下Huber 估计器和 WLS 的差别不到 5%攻击者只要改两个关键量测就能把状态拉偏。后来补了 PMU 的量测冗余度提到 2.1同样的攻击场景下 Huber 估计器的状态偏差从 12% 降到 3% 以内。所以做鲁棒估计之前先算冗余度低于 1.8 的话优先补量测而不是调算法。3. 用 Python 跑通一个最小鲁棒状态估计器3.1 构造 IEEE 14 节点算例与量测向量先搭一个能复现的算例。用 pandapower 建 IEEE 14 节点模型生成潮流真值再按真值加高斯噪声造量测最后注入虚假数据。这一步的关键是攻击向量要满足 FDIA 的构造条件即攻击后的量测残差与攻击前几乎一致否则随便加个噪声都能被检测出来测不出鲁棒估计器的真实能力。import numpy as np import pandapower as pp import pandapower.networks as pn # 建 IEEE 14 节点模型并跑潮流 net pn.case14() pp.runpp(net) # 提取真值节点电压幅值、相角支路功率 V_true net.res_bus.vm_pu.values theta_true np.deg2rad(net.res_bus.va_degree.values) P_branch_true net.res_line.p_from_mw.values Q_branch_true net.res_line.q_from_mw.values # 组装量测向量 z h(x) e这里简化为直接用量测函数 # 实际工程中 h(x) 是非线性潮流方程这里用真值加噪声模拟 np.random.seed(42) sigma 0.01 # 量测噪声标准差 z_voltage V_true np.random.normal(0, sigma, len(V_true)) z_power P_branch_true np.random.normal(0, sigma * 100, len(P_branch_true)) # 构造 FDIA 攻击向量攻击者篡改 3 号和 8 号节点的注入功率量测 # 攻击量 a 满足 a H * cc 是状态偏移向量H 是量测雅可比矩阵 # 这里简化处理直接在量测上加一个与拓扑相关的偏移 attack_idx [2, 7] # 对应节点 3 和 8 z_power_attacked z_power.copy() z_power_attacked[attack_idx] np.array([15.0, -12.0]) # 注入虚假功率偏移 print(f攻击前量测均值: {z_power.mean():.2f}) print(f攻击后量测均值: {z_power_attacked.mean():.2f})这段代码做了三件事跑潮流拿真值、加噪声造量测、在指定量测上注入偏移。参数sigma控制噪声水平实际工程中功率量测的噪声标准差通常在 1%~2% 额定值这里用sigma * 100是因为功率基准是 100MW。攻击偏移量 15MW 和 -12MW 是随手设的真实攻击者会按a Hc构造让残差不变但这里为了演示鲁棒估计器的压制效果直接用固定偏移就够了。3.2 Huber 估计器的迭代重加权实现Huber 估计器的核心是迭代重加权最小二乘IRLS。每一轮用当前残差算权重残差大的量测权重低然后解一次 WLS更新状态再算残差直到收敛。下面是一个简化版实现状态变量只取电压幅值和相角量测函数用线性化近似。def huber_weight(residual, delta1.5): Huber 权重函数残差小于 delta 时权重为 1超过时按 delta/|r| 衰减 abs_r np.abs(residual) weights np.ones_like(abs_r) mask abs_r delta weights[mask] delta / abs_r[mask] return weights def robust_state_estimation(z, H, x0, max_iter20, tol1e-6, delta1.5): z: 量测向量 (m,) H: 量测雅可比矩阵 (m, n) x0: 状态初值 (n,) delta: Huber 阈值通常取 1.5 倍量测标准差 x x0.copy() for it in range(max_iter): # 计算残差 r z - H x # 算 Huber 权重 w huber_weight(r, delta) # 加权最小二乘解x (H^T W H)^(-1) H^T W z W np.diag(w) HtWH H.T W H HtWz H.T W z x_new np.linalg.solve(HtWH, HtWz) # 收敛判断 if np.linalg.norm(x_new - x) tol: print(f收敛于第 {it1} 次迭代) break x x_new return x, w # 构造简化的量测雅可比矩阵实际应用需按潮流方程求偏导 # 这里用随机矩阵模拟仅演示算法流程 m, n 20, 10 H np.random.randn(m, n) x_true np.random.randn(n) z_clean H x_true np.random.normal(0, 0.01, m) z_attack z_clean.copy() z_attack[2] 0.5 # 注入攻击 z_attack[7] - 0.4 x0 np.zeros(n) x_est_clean, w_clean robust_state_estimation(z_clean, H, x0) x_est_attack, w_attack robust_state_estimation(z_attack, H, x0) print(f干净数据状态误差: {np.linalg.norm(x_est_clean - x_true):.4f}) print(f攻击数据状态误差: {np.linalg.norm(x_est_attack - x_true):.4f}) print(f被攻击量测的权重: {w_attack[2]:.3f}, {w_attack[7]:.3f})这段代码里huber_weight是权重函数delta是阈值取 1.5 倍量测标准差是经验值。robust_state_estimation做 IRLS 迭代每次用当前残差更新权重再解加权最小二乘。关键参数max_iter控制最大迭代次数tol是收敛容差。运行后你会看到干净数据下状态误差很小攻击数据下误差被压制而且被攻击量测的权重明显低于 1。这就是鲁棒估计器在起作用——它没有去识别哪个量测被攻击而是通过降权让攻击量测无法主导解。3.3 用残差协方差做攻击检测的辅助判据鲁棒估计器本身不输出有没有攻击的结论它只是让估计结果更稳。如果你需要报警还得加一个检测环节。常用的是归一化残差检验算每个量测的残差除以其标准差超过阈值就报警。但 FDIA 的残差可能很小所以更可靠的是用鲁棒估计器的权重分布——如果大量量测权重同时下降说明系统里存在系统性偏差而不是单个坏数据。def attack_detection(w, threshold0.5, ratio0.3): w: 鲁棒估计器输出的权重向量 threshold: 权重低于此值视为可疑 ratio: 可疑量测比例超过此值触发报警 suspicious np.sum(w threshold) suspicious_ratio suspicious / len(w) if suspicious_ratio ratio: return True, suspicious_ratio return False, suspicious_ratio # 用上面的权重做检测 is_attack_clean, ratio_clean attack_detection(w_clean) is_attack, ratio_attack attack_detection(w_attack) print(f干净数据可疑比例: {ratio_clean:.2%}, 报警: {is_attack_clean}) print(f攻击数据可疑比例: {ratio_attack:.2%}, 报警: {is_attack})这个检测逻辑很简单统计权重低于 0.5 的量测比例超过 30% 就报警。参数threshold和ratio需要根据实际系统的量测冗余度和噪声水平调。冗余度高的系统可以放宽ratio因为正常量测多少数被降权不影响比例。冗余度低的系统要收紧否则容易漏报。4. 避坑与排查鲁棒状态估计器落地时的五个血泪教训4.1 现象估计结果震荡不收敛迭代 50 次还在跳原因Huber 阈值delta设得太小正常量测也被降权权重矩阵每轮剧烈变化IRLS 在解附近来回震荡。或者量测雅可比矩阵H的条件数太大加权后更病态。解决先把delta调到 2.0~2.5 倍量测标准差观察收敛曲线。如果还震荡检查H矩阵的条件数超过 1e6 的话需要做量测筛选或加正则化项。我一般会在HtWH上加一个小的对角项1e-6 * I相当于岭回归能显著改善数值稳定性。4.2 现象攻击量测的权重没降下来估计结果还是被带偏原因攻击者构造的虚假数据与正常量测的残差分布很接近Huber 权重函数在阈值附近区分度不够。或者攻击量测的数量超过了鲁棒估计器的崩溃点breakdown pointHuber 的崩溃点约 50%但实际有效范围通常只有 30%。解决换 Tukey 双权函数它的红降特性对接近阈值的残差更敏感。或者引入投影统计量做预筛把可疑量测在迭代前就标记出来给它们更低的初始权重。如果攻击量测比例确实超过 30%单靠鲁棒估计器不够需要结合 PMU 的动态量测做交叉验证。4.3 现象投影统计量计算太慢实时性达不到要求原因投影统计量要对每个量测计算所有投影方向的中位数和 MAD复杂度 O(m²)m 是量测数。区域电网 m 可能上千单次计算就超过状态估计本身的时间。解决只对残差最大的前 20% 量测做投影统计其余量测直接给正常权重。或者用随机投影代替全方向投影随机选 50~100 个方向精度损失很小但速度提升一个数量级。另一个做法是把投影统计量放在状态估计之前做一次后续迭代不再重复计算。4.4 现象量测冗余度不足时鲁棒估计器和 WLS 结果几乎一样原因冗余度低于 1.5 时坏数据检测的自由度不够鲁棒估计器的权重调整空间被压缩。攻击者只要改少数关键量测就能同时骗过 WLS 和鲁棒估计器。解决优先补量测尤其是 PMU 的电压相角量测它对状态估计的可观测性贡献最大。如果补不了量测退而求其次用历史数据做时序一致性检验——攻击者可以改单点量测但很难同时改多个时间断面的量测而保持时序连贯。把时序残差也纳入权重计算能部分弥补冗余度不足。4.5 现象攻击检测误报率高正常操作也被报警原因检测阈值ratio设得太低或者系统本身存在量测偏差比如 CT 慢漂移导致正常量测的权重也偏低。另外如果系统里有大量零注入节点这些节点的量测权重天然不稳定容易触发误报。解决先做一轮无攻击场景的基线测试统计正常情况下的可疑量测比例把ratio设成基线的 2~3 倍。对零注入节点单独处理不纳入可疑比例统计。如果 CT 漂移是已知问题在状态估计之前先做量测校准别让鲁棒估计器去扛这个锅。5. 从离线验证到在线部署一个可复用的验证套路鲁棒状态估计器写完只是第一步怎么证明它在真实攻击下有效才是决定要不要投入的关键。我一般会走三步验证离线注入测试、半实物仿真、现场试运行。离线测试用历史量测数据人为注入不同比例的 FDIA看状态偏差和检测率。半实物仿真用 RTDS 或者 RT-LAB 接真实 PMU验证通信延迟和量测丢包对鲁棒估计器的影响。现场试运行先旁路运行不接入闭环控制只记录估计结果和报警日志跑两周再评估。下面是一个离线验证的脚本框架用蒙特卡洛跑 100 次不同攻击场景统计状态误差和检测率。def monte_carlo_validation(n_trials100, attack_ratio0.2): n_trials: 蒙特卡洛次数 attack_ratio: 被攻击量测的比例 errors_wls [] errors_huber [] detection_rates [] for trial in range(n_trials): # 每次重新生成量测和攻击 m, n 30, 12 H np.random.randn(m, n) x_true np.random.randn(n) z H x_true np.random.normal(0, 0.01, m) # 随机选 attack_ratio 比例的量测注入攻击 n_attack int(m * attack_ratio) attack_idx np.random.choice(m, n_attack, replaceFalse) z_attack z.copy() z_attack[attack_idx] np.random.normal(0, 0.5, n_attack) # WLS 估计 x_wls np.linalg.lstsq(H, z_attack, rcondNone)[0] errors_wls.append(np.linalg.norm(x_wls - x_true)) # Huber 估计 x_huber, w robust_state_estimation(z_attack, H, np.zeros(n)) errors_huber.append(np.linalg.norm(x_huber - x_true)) # 检测 is_attack, _ attack_detection(w) detection_rates.append(1 if is_attack else 0) print(fWLS 平均状态误差: {np.mean(errors_wls):.4f}) print(fHuber 平均状态误差: {np.mean(errors_huber):.4f}) print(f攻击检测率: {np.mean(detection_rates):.2%}) print(f误差降低幅度: {(1 - np.mean(errors_huber)/np.mean(errors_wls)):.2%}) monte_carlo_validation(n_trials100, attack_ratio0.2)这个脚本跑 100 次每次随机选 20% 的量测注入攻击对比 WLS 和 Huber 的状态误差。参数attack_ratio可以调从 0.1 到 0.4 各跑一遍看鲁棒估计器的误差降低幅度怎么变化。如果attack_ratio超过 0.3 后误差降低幅度骤降说明这个配置下的崩溃点到了需要补量测或者换更强的鲁棒估计器。验证通过后在线部署还有几个工程细节要注意。第一状态估计的周期通常是 5~15 秒鲁棒估计器的迭代次数要控制在这个时间窗内max_iter别超过 10。第二权重矩阵W的存储和计算要优化用稀疏矩阵别用稠密np.diag。第三报警日志要记录每次迭代的权重分布方便事后回溯——攻击者可能慢慢调大量测偏移单次看不出来但权重分布的趋势会暴露问题。我自己的习惯是每次现场试运行前先用历史数据跑一遍离线验证把delta、ratio、max_iter这三个参数记在配置文件里别硬编码在代码里。现场环境一变量测噪声水平可能差一倍参数不调的话鲁棒估计器要么不收敛要么误报率飙升。这套东西没有一劳永逸的参数只有不断根据现场数据微调的习惯。希望帮到你。本文还有配套的精品资源点击获取
