简介面向电力系统状态估计与网络攻击防御研究者的MATLAB源码包聚焦基于鲁棒广义极大似然GM估计器的虚假数据注入攻击防御方法。方法融合投影统计与Givens旋转可同时抵御坏数据、坏杠杆点、坏零注入及相关网络攻击并在高斯或厚尾非高斯噪声下保持较高统计效率。压缩包共13个文件包含10个m源程序核心估计器、IEEE标准算例数据、导纳矩阵计算、量测转换与GM/WLS对比测试模块、1份PDF与1份DOCX说明文档、1个txt说明压缩包整体仅159KB轻量易部署。已有1104人浏览学习。通过运行测试脚本可直观复现GM估计器在受攻击量测下的鲁棒性并支持变压器抽头位置联合估计适合作为科研对比基准、课程设计拓展或工程验证参考。1. 调度大屏上的量测曲线在悄悄漂移鲁棒状态估计器到底挡的是什么调度大屏上的量测曲线在悄悄漂移EMS的坏数据检测却毫无反应——这不是量测噪声而是虚假数据注入攻击FDIA在绕开残差检验。攻击者只要能从公开的电网参数里算出量测雅可比矩阵就能设计一个攻击向量让基于加权最小二乘WLS的状态估计在新的错误运行点上正常工作。基于鲁棒电力系统状态估计器的防御方法核心是把目标函数从二次型换成对极端残差不敏感的形式让注入数据在目标函数值和权重分布上暴露出来。这篇笔记面向能跑MATLAB或Python、手上有IEEE节点数据的工程师和学生把原理、可复现的攻防实验、关键参数和落地踩坑一次性讲清楚。2. 从WLS到鲁棒估计FDIA攻击面与目标函数的选型逻辑2.1 传统WLS的最优性依赖高斯假设电力系统状态估计的量测集合来自SCADA和PMU节点注入有功无功、支路潮流、母线电压幅值。模型统一写成 z h(x) ex 是状态量通常是各节点电压幅值和相角h(x) 由潮流方程决定e 是量测误差工程上近似为零均值高斯噪声。加权最小二乘WLS的做法是把加权残差平方和 J(x) 最小化权重矩阵 W 取各量测方差的倒数对角阵。一阶最优性条件是 H^T W [z - h(x)] 0其中 H 就是 h(x) 对 x 的雅可比矩阵实际求解用高斯-牛顿法迭代每一步解一个 (H^T W H) Δx H^T W r 的线性方程组。WLS 能在电力系统里统治几十年是因为量测噪声近似高斯时最小二乘在无偏和最小方差意义上是最优的。坏数据检测也建立在这个分布假设上状态估计算完以后目标函数 J(x̂) 在正常工况下服从自由度 m-n 的卡方分布取卡方分布 95% 或 99% 分位数当阈值超过就怀疑有坏数据再配合最大归一化残差逐点定位。具体来说残差 r z - h(x̂) 的协方差矩阵是 S W^{-1} - H (H^T W H)^{-1} H^T归一化残差 r_i^N |r_i| / sqrt(S_ii)通常超过 3.0 就标记为可疑。这套检测机制的命门在于一个硬假设异常必须表现为残差变大。坏数据检测器看到的永远只是残差大不大它不关心状态量本身合不合理。后面的内容全部围绕一个反命题展开——攻击者能不能让残差完全不变同时把状态估计结果推到任意想要的错误位置。答案是能而且数学上非常干净。2.2 FDIA如何把攻击向量藏在雅可比矩阵的列空间虚假数据注入攻击从 2011 年前后被系统研究核心思想是生成一套与残差检测完全自洽的量测篡改。攻击者选一个任意非零向量 c令攻击向量 a H c然后把量测改成 z_a z a。在直流线性化模型 z H x e 下这个结论一目了然x̂ (H^T W H)^{-1} H^T W zx̂_a (H^T W H)^{-1} H^T W (z Hc) x̂ cr_a z_a - H x̂_a z Hc - H(x̂ c) z - H x̂ r也就是说攻击后的状态估计结果变成了 x̂ c而残差和攻击前一模一样。卡方统计量不变最大归一化残差不变任何基于残差的统计检验全部失效。交流模型下严格说攻击向量需要按攻击后的运行点重新计算但只要攻击的 c 取小一些在局部运行点用 H(x̂) c 构造注入一阶近似误差就小到检测器完全无感。按攻击者掌握的信息量FDIA 可以分三档理解完全知晓型掌握完整 H 矩阵能构造严格无残差攻击这是防御研究默认的最坏假设部分知晓型只知道部分线路参数构造的攻击向量不完全落在列空间里残差会有微量变化小幅度注入仍可能被阈值放过零知识盲攻击只能拿公开数据统计推断 H效果取决于信息暴露程度。做防御评估时应该固定用第一档能挡住完全知晓型攻击低阶攻击自然失效。还有一个工程现实攻击者希望篡改的量测数量尽可能少以减少被物理机制发现的风险。但稀疏篡改和a 严格等于 Hc存在矛盾现实中攻击者常通过求解稀疏优化问题来近似结果就是实际攻击的残差并非严格为零。这个非零量给了检测器操作空间也是后面鲁棒估计器的阈值参数能起作用的前提。2.3 M估计、LAV与GM估计换目标函数为什么能破局既然二次型目标函数是根子上的弱点自然想到换一个对异常残差不敏感的目标函数。M估计把目标函数写成 min Σ ρ(r_i / σ)ρ 是增长比平方更慢的函数。最常用的是 Huber 函数ρ(r) r²/2, |r| ≤ δρ(r) δ|r| - δ²/2, |r| δ对应的等价权函数是 w(r) ρ(r) / r落在 δ 之内权重为 1超出后权重按 δ/|r| 衰减。这样一两个被严重篡改的量测无论残差多大对目标函数的贡献从平方级降到线性级状态估计解就不会被整体拖跑。求解用迭代重加权最小二乘IRLS每一步把当前权重代入标准 WLS 框架实现成本低和现有 WLS 代码天然兼容。比 M 估计更刚的是最小绝对误差LAV估计目标函数取 min Σ |r_i|。LAV 的崩溃点理论上是 50%比 M 估计更能扛高比例坏数据但目标函数不可微线性化后要写成线性规划求解变量维度随系统规模膨胀大电网里计算成本明显更高。还有一种广义 M 估计GM 估计在 M 估计基础上额外计算每个量测的杠杆值帽子矩阵对角元大的量测被额外降权专门对付既带异常残差又处在高杠杆位置的量测点。三种估计器的选型逻辑我按落地场景给三条想快速替换现有 WLS 代码、迭代框架越小改动越好选 Huber M 估计仿真里攻击比例很高或者量测冗余度不足选 LAV系统存在天然杠杆点、量测冗余度分布严重不均选 GM 估计。实际电网里高连接度母线的注入量测就是天然杠杆点攻击者优先盯这些位置所以 FDIA 防御场景下 GM 估计最贴合工程现实这一点在避坑章还会展开。3. 最小可复现实验搭一套FDIA攻防测试平台3.1 平台选型MATPOWER vs pandapower先把用什么跑通这件事定下来。MATLAB 下的 MATPOWER 内置了交流状态估计和坏数据检测开箱即用但要构造 FDIA 攻击向量、提取雅可比矩阵、自定义目标函数这些攻防实验的核心环节反而要自己写MATLAB 授权和服务器批量跑批也是个现实约束。我更推荐 Python 路线pandapower 负责加载 IEEE 节点系统和生成量测numpy/scipy 承担矩阵运算与优化整个流程在 Linux 服务器上纯命令行跑批非常顺。两者的取舍可以看下面这张表。对比维度MATLAB MATPOWERPython pandapower状态估计能力内置 run_se开箱即用需要自己装配迭代或装额外模块自定义雅可比与攻击构造需自行写脚本矩阵访问方式偏 MATLAB 习惯数组操作直接numpy 生态统一批量实验与服务器部署依赖 MATLAB 运行时多节点跑批麻烦纯命令行无授权负担数据可视化图窗交互方便matplotlib 多两行代码但可脚本化我的习惯是用 pandapower 加载网络数据但 FDIA 构造和鲁棒估计器全部用 numpy 手写不依赖工具包封装好的状态估计接口。这样每一步矩阵运算都是显式的出了数值问题好定位。你手上那个 zip 包无论里面是 MATLAB 还是 Python 代码落地最快的方式都是先跑通一个最小直流线性化实验理解残差为什么不变再往交流牛顿法上靠。3.2 构造FDIA攻击向量一个可以直接跑的Python片段下面这段代码把 FDIA 攻防闭环浓缩到最小规模。我用一个随机生成的对称正定矩阵代替节点电纳矩阵去掉平衡节点后得到雅可比 H量测真值由 H θ 生成再叠高斯噪声。攻击者构造 c 并计算 a H c把攻击加进量测然后分别跑攻击前后的 WLS输出最大归一化残差。import numpy as np rng np.random.default_rng(42) n_state 14 # 状态量数量相角对应IEEE 14节点规模 sigma 0.02 # 量测噪声标准差单位弧度 # 构造对称正定矩阵模拟电纳阵保证去掉平衡节点后 H 满列秩 A rng.normal(size(n_state, n_state)) B 5.0 * A A.T 6.0 * np.eye(n_state) theta_true np.array([0.00, 0.05, -0.03, 0.02, 0.01, -0.04, 0.03, -0.02, 0.06, -0.01, 0.00, 0.02, -0.05, 0.04]) # 单位弧度 H B[1:, 1:] # 平衡节点相角固定为0去掉第0行/列 theta_est theta_true[1:] # 生成正常量测H theta 高斯噪声 z H theta_est rng.normal(0.0, sigma, sizeH.shape[0]) def wls_solve(H, z): G H.T H return np.linalg.solve(G, H.T z) x_wls wls_solve(H, z) r_norm z - H x_wls print(攻击前最大归一化残差: %.4f % np.max(np.abs(r_norm) / sigma)) # 构造FDIAc是攻击者希望注入的状态偏移攻击向量 a H c c np.full(H.shape[1], 0.03) # 每个状态量偏移0.03 rad可改成稀疏模式 a H c z_att z a x_att wls_solve(H, z_att) r_att z_att - H x_att print(攻击后最大归一化残差: %.4f % np.max(np.abs(r_att) / sigma))逻辑说明wls_solve 对直流模型求闭式解不需要迭代。先算攻击前的最大归一化残差正常噪声下这个值在 23 之间攻击后残差和攻击前几乎相等这正是 FDIA 绕过卡方检验的直接证据。c 取了统一偏移 0.03 弧度真实攻击会设计成只篡改部分量测的稀疏形式但数学本质相同。换到真实 IEEE 14 节点时只需要把 B 换成由线路电抗组成的节点电纳矩阵其他代码不用动。参数说明sigma 取 0.02 弧度对应工程上量测误差的典型水平取太小会让攻击检测率虚高取太大会让正常噪声本身就掩盖攻击。随机种子固定为 42 保证可复现。这里的 H 是 13×13 的稠密矩阵真实电网的 H 是稀疏的但 FDIA 原理不依赖稀疏性用稠密矩阵做演示更直观。注意直流近似模型下攻击向量严格满足 a Hc交流模型只有在攻击偏移较小时近似成立。调试攻防流程时永远先用直流模型跑通再切交流牛顿法。3.3 把实验跑在Linux服务器上解压、目录与运行拿到 zip 包的第一件事不是看代码而是确认环境。Linux 下解压一行命令unzip 包名.zip。我习惯先执行 unzip -t 做完整性检查遇到 missing zip entry 或者 end-of-central-directory signature not found 这类报错十有八九是传输丢字节或下载中断重新传一遍比抠参数高效。伪加密的问题也常在这个时候冒出来——解压提示输入密码其实文件根本没加密用 7z 等工具可以直接绕过伪加密标志位解出来。# 解压并确认目录结构 unzip -t 鲁棒电力系统状态估计器FDIA防御.zip unzip -o 鲁棒电力系统状态估计器FDIA防御.zip -d fdia_attack cd fdia_attack python -m venv venv source venv/bin/activate pip install numpy scipy pandapower matplotlib# 最小化验证跑通线性化攻防演示 python scripts/fdia_demo.py逻辑说明第一个代码块做了三件事——完整性检查、解压到指定目录、创建虚拟环境并安装依赖。把依赖隔离到 venv 里能避免污染系统 Python。第二个代码块运行演示脚本输出两条归一化残差。如果两条值几乎一样说明 FDIA 构造成功如果攻击后残差明显变大说明 H 的列空间没对上最常见的错误是平衡节点没剔除干净导致 H 不满列秩。参数说明unzip -o 表示覆盖同名旧文件-d 指定输出目录。python -m venv 创建当前目录下的虚拟环境source 激活后 pip 安装的包都落在 venv 里不碰系统环境。numpy 和 scipy 是矩阵运算基础pandapower 只在需要加载真实 IEEE 节点数据时用matplotlib 画攻防对比曲线。离线内网环境就把依赖包下载成 whl 文件再离线安装本质和普通 pip 流程一样。4. 鲁棒估计器的关键参数残差尺度、Huber阈值与迭代收敛4.1 残差尺度估计1.4826这个系数是谁定的M 估计目标函数里的残差必须先除以尺度因子 σ 才能无量纲化σ 的选法直接决定权重函数的形状。σ 偏大估计器对注入攻击的灵敏度下降σ 偏小正常噪声被误判成坏数据估计效率受损。最稳的做法不靠人为设定而是用中位数绝对偏差MAD从残差样本里直接估计σ 1.4826 × median(|r_i - median(r)|)系数 1.4826 的来历是高斯分布下 MAD 与标准差的比值约为 0.6745取倒数就是 1.4826。这样估计出的 σ 在纯高斯噪声下逼近真实标准差而且对个别离群值不敏感正好配合鲁棒估计器的迭代过程。实现上每轮迭代都要重新用当前残差估计 σ而不是只在初始时算一次。需要注意一个边界当被污染的量测比例超过一半时MAD 本身也会被拉偏。这时候两个处理办法一是改用更鲁棒的尺度估计量比如分位数间距二是在迭代中把 σ 固定在攻击前标称量测噪声水平上。后者在防御场景里更实用因为调度系统对量测精度的标称值是有台账的拿台账值兜底比纯数据驱动更可靠。4.2 Huber阈值δ1.345是个好起点但不是定死的神话Huber 函数的分段点 δ 控制多大残差算异常。δ 越大估计器越接近 WLS抗攻击能力弱但正常段效率高δ 越小抗攻击能力强但正常噪声下的方差变大。统计文献里经典建议是 δ 1.345σ这个值使估计器在纯高斯噪声下相对 WLS 的效率保持在 95% 左右。但防御 FDIA 和做理论统计不是一回事我习惯把 δ 放在 1.0σ 到 2.0σ 之间做敏感度扫描再根据攻击强度和误报率的平衡定值。def huber_weights(r, sigma, delta): Huber等价权残差超过 delta*sigma 的量测线性降权 rs np.abs(r) / sigma w np.ones_like(rs) mask rs delta w[mask] delta / rs[mask] return w def robust_wls(H, z, sigma, delta1.345, max_iter30, tol1e-6): IRLS求解Huber M估计 x np.linalg.lstsq(H, z, rcondNone)[0] # WLS解作初值 for _ in range(max_iter): r z - H x w huber_weights(r, sigma, delta) W np.diag(w) x_new np.linalg.solve(H.T W H, H.T W z) if np.max(np.abs(x_new - x)) tol: return x_new, w, r x x_new return x, w, r逻辑说明huber_weights 实现等价权函数残差超过 δσ 的量测权重从 1 掉到 δ/|rs|残差越大权重越小但不会归零这保证被攻击量测仍保留少量信息而不是被彻底丢弃。robust_wls 是完整的 IRLS 迭代WLS 解做初值每轮重算残差、更新权重、解加权正规方程。初值用 np.linalg.lstsq 而不是全零是为了避免迭代落到差的局部解。参数说明delta 需要和 sigma 配套理解代码里实际阈值是 delta 乘以 sigma。max_iter 给 30 次中小规模系统通常 510 次收敛。tol 取 1e-6按状态量最大增量判停。这个实现对条件数较差的 H 会暴露数值问题条件数超过 1e8 时先做列归一化再用。4.3 收敛判据看状态量还是看权重变化IRLS 迭代停在哪直接反映鲁棒估计器是否真正稳定。我同时看两个信号状态量的最大增量小于 1e-4同时所有权重中最大变化量小于 1e-3两者都满足才认为收敛。只盯状态量容易在权重还在振荡时提前停只盯权重又可能在状态量未收敛时误判。实际仿真遇到迭代发散优先查 H 矩阵条件数条件数差到 1e12 时调参没用要换量测组合或者改求解方法。参数推荐起点调整方向备注残差尺度 σ1.4826 × MAD(r)污染比例高时固定为标称噪声高斯下逼近真实标准差Huber δ1.345σ抗攻击优先取小效率优先取大小于 1.0σ 会误伤正常量测迭代上限30 次条件数差时加到 50一般 510 次收敛状态量收敛阈值1e-4更高精度取 1e-6与权重判据同时使用权重变化阈值1e-3防止权重未稳定时早停单独使用会误判这张参数表是多次攻防仿真沉淀出来的起点不是标准答案。电网规模、量测冗余度、攻击比例都会影响最优取值第 6 章会给一套完整的扫描验证方法避免单点调参带来的侥幸结果。5. 避坑记录鲁棒状态估计器落地时最容易翻车的5件事5.1 稀疏攻击把残差检验全骗过去了现象用 FDIA 攻击向量替换部分量测后最大归一化残差和卡方统计量与攻击前几乎一致传统坏数据检测完全失效WLS 解被推动了几十兆瓦但系统毫无告警。原因攻击向量 a Hc 落在雅可比矩阵的列空间里攻击后残差在数学上不变。任何只看残差表面大小的检测器都看不到它这不是检测阈值设得松而是检测量本身被攻击者绕过了。解决把状态估计器换成 Huber M 估计让鲁棒估计器和传统残差检测双通道运行。鲁棒估计器对这类攻击的响应是权重显著下降、目标函数值异常抬升用这两个信号做判据就能把 FDIA 和普通坏数据区分开。实际项目里我在鲁棒估计器后端串了一个告警逻辑权重跌破 0.3 的量测数量超过预设比例就触发报警效果比单纯盯目标函数更直观。5.2 电压相角的量纲不统一残差对比是假象现象攻击前后残差都很大甚至攻击后残差比攻击前还小看起来 FDIA 无效实际是构造攻击时把角度和幅值混在同一个数量级里。原因交流状态估计里相角单位是弧度、电压幅值单位是标幺值两者量级差几十倍。H 矩阵不同列的数值差异巨大正规方程条件数恶化攻击向量也失真残差对比完全不可信。解决所有量测先做无量纲化角度用弧度、幅值用标幺值残差按各自量测精度归一化后再进权重函数。直流近似模型天然没有这个问题所以调试阶段用直流模型是对的切交流模型前先把量纲检查清单走一遍能省下大量排查时间。5.3 IRLS迭代不收敛目标函数值在震荡现象Huber 估计器在 IEEE 30 节点系统上迭代十几轮不收敛目标函数值减不下去状态量在两个值之间来回跳。原因权重矩阵每一步完全重算当坏数据比例接近 50% 时等价权更新步长过大迭代进入极限环H 矩阵病态时权重变化被数值误差进一步放大。解决给权重更新加阻尼状态量更新时乘一个阻尼因子 λ0.5 再与旧值合并或者限制单次权重变化幅度同时检查 H 条件数超过 1e8 做列归一化。血泪经验是先处理条件数再加阻尼顺序反了容易把数值问题误判成算法不收敛。5.4 M估计对杠杆点攻击不鲁棒估计值照样被拉偏现象攻击加在某个高连接度母线的注入量测上Huber 估计器的残差统计量没有明显异常但状态估计结果仍然被拉偏权重函数完全没触发。原因M 估计只按残差大小降权而杠杆点位于 H 阵的行空间边缘正常工况下杠杆点的残差本来就不大。攻击者只要沿着杠杆方向构造状态偏移残差照样很小纯残差型 M 估计对它没有反应。解决升级成 GM 估计在权重计算里乘一个杠杆惩罚因子帽子矩阵对角元 h_ii 大的量测被自动降权。实现上有现成的公式可抄w_lev sqrt(1 - h_ii)和 Huber 等价权相乘作为最终权重。这个坑在纯理论仿真里容易被忽略实际电网量测分布不均时很常见。5.5 攻击强度设置不合理导致防御指标虚高或虚低现象仿真报告里检测率 99%换一个攻击强度或随机种子就掉到 60%指标完全不可复现结论没法让人信服。原因攻击强度太大注入量远超正常量测波动任何鲁棒估计器都能抓出来检测率虚高攻击强度太小攻击淹没在量测噪声里再好的算法也白搭。单点调参测试的本质是在噪音里碰运气。解决把攻击强度定义为量测噪声标准差的倍数 α核心指标做成 α 从 0.5 到 5 的扫描曲线同时固定多个随机种子取平均。这样得到的检测率曲线有实际参考价值也避免了对某个攻击参数的偶然依赖。这条算是我在攻防实验里踩得最深的一个坑。6. 进阶验证技巧用受控强度攻击评估防御有效性6.1 用攻击强度扫描曲线替代单点测试鲁棒估计器到底扛得住多强的注入单点测试说明不了问题。常见做法是把攻击向量写成 a α × (H c)α 是攻击强度倍数c 固定为单位范数的状态偏移然后让 α 从 0.5 扫到 5.0。对应 0.02 弧度的量测噪声α1 表示注入偏移和噪声同量级α5 是显著恶意攻击。对每个 α 跑多个随机种子把 WLS 和鲁棒估计器的状态偏差取平均画两条曲线。鲁棒估计器曲线在 α 增大时增长明显平缓说明防御有效两条曲线几乎重合说明 δ 或 σ 配置有问题回到第 4 章调参数。6.2 评估指标检测率、误报率与估计偏差攻防实验结论要落到三个指标上检测率 DR、误报率 FAR、状态估计偏差。DR 和 FAR 衡量检测本身估计偏差衡量鲁棒估计器在攻击下的性能保持三者放同一张表里看才有意义。评估项计算方法参考目标检测率 DR报警数 / 攻击数多随机种子平均α≥2 时趋于 95% 以上误报率 FAR正常量测误警数 / 正常实验数控制在 5% 以下状态估计偏差‖x̂ - x_true‖₂ 的均值与标准差攻击下不显著恶化权重量变化比例攻击前后权重明显下降的量测占比与攻击稀疏模式匹配我早期吃过大亏只在一个攻击强度下测出 99% 的检测率就把结论写进报告换到 α1 左右立刻翻车。后来固定了扫描曲线加多随机种子的验证流程指标才算稳定下来。现在做任何对比实验先跑一遍扫描曲线再看单点指标。6.3 一个值得保留的调试习惯调试鲁棒估计器时我习惯把每次迭代的权重向量打出来而不是只盯最终检测率。权重按量测索引画出来被攻击的量测会明显变成深色如果攻击后权重没有明显变化说明阈值 δ 太大或 σ 估计失真算法根本没进入鲁棒模式。这个可视化习惯帮我排掉过好几个隐蔽的数值错误比如 H 矩阵代入时算错、σ 固定值取错数量级。你搭实验平台时不妨把每轮的 w 都存下来跑通了再删调试代码。希望帮到你。本文还有配套的精品资源点击获取
