1. 项目概述：为什么电力系统状态估计要引入鲁棒滤波
做电力系统动态状态估计的人，十有八九都跟扩展卡尔曼滤波器（EKF）打过交道。但真到了实际工程里，EKF的表现往往没有仿真论文里那么漂亮——最大的痛点就是对异常值（outlier）太敏感。量测里只要混进去一个坏数据，比如PMU通信丢包产生的跳变值、互感器饱和瞬间的毛刺，滤波器的状态估计结果立刻就被带偏，甚至直接发散。
这个项目做的就是把经典的EKF往两个方向改造：一是引入Huber鲁棒化代价函数，降低异常量测对状态更新的影响权重；二是在每个采样时刻做迭代更新，缓解强非线性系统里一阶线性化带来的误差累积。两个方向结合在一起，就是标题里说的“鲁棒迭代扩展卡尔曼滤波器”（Robust Iterative Extended Kalman Filter，RIEKF）。
这套滤波器适合谁来用？我认为三类人最需要：第一类是研究电力系统动态状态估计的硕博生，需要跟EKF、UKF、鲁棒H∞滤波做对比实验的；第二类是做在线状态监测的工程技术人员，面对的量测数据质量参差不齐，需要抗差能力的估计器；第三类是刚入门滤波理论、想搞明白“鲁棒化到底改在哪里”的学习者。项目用Matlab实现，代码结构清晰，改参数方便，跑通一遍基本就能把RIEKF的脉络摸清楚。
整个实现里最核心的思路，是把鲁棒化和迭代化分别作用在滤波器的不同环节：Huber方法修正的是“量测更新”里的增益计算，迭代机制修正的是“状态更新”里的线性化精度。两个改动各有各的数学依据，也各有各的适用边界，下面我把设计思路和代码实现逐层拆开讲。
2. 核心原理拆解：EKF为什么需要“迭代+鲁棒”双重改造
2.1 EKF的线性化短板与迭代补偿思路
EKF处理非线性系统的基本思路是围绕当前估计点做一阶泰勒展开。每来一个新量测，就先算雅可比矩阵，把非线性方程近似成线性方程，再套标准卡尔曼滤波更新公式。这个思路在系统非线性较弱时没问题，但电力系统的动态模型往往非线性程度不低——发电机功角方程里的正弦项、励磁系统里的限幅环节，一阶近似误差就会随着时间步长累积。
迭代扩展卡尔曼滤波器（IEKF）的做法，是把“线性化—更新—再线性化”这个循环在一个采样周期内重复执行若干次。每更新一次状态，就用新状态重新计算雅可比矩阵，再次执行量测更新。这个过程本质上是在求解一个高斯-牛顿类型的优化问题，迭代次数足够多时，结果会逼近最大后验估计（MAP）的数值解。
我做个生活化类比：EKF相当于你问路时只听了第一个路人的指路，然后闷头往前走；IEKF则是每走完一小段就重新问一次路，路线自然更准。代价就是计算量翻倍，但换来的精度提升在强非线性系统里非常可观。
2.2 Huber鲁棒化：给异常值“降权”的数学方案
Huber鲁棒化的核心思想并不复杂——当残差（实际量测与预测量测之差）过大时，不要继续按平方项处理它。标准EKF的代价函数是残差的二次型，大残差会被平方放大，导致滤波器疯狂追赶异常值。Huber函数则把代价函数设计成分段形式：
- 当残差小于阈值时，按二次型处理，保持最优性；
- 当残差大于阈值时，按一次型处理，把异常值的影响限制在线性增长而不是二次爆炸。
阈值通常选取为1.345倍的残差标准差，这个数值在Huber的原始文献里给出了理论依据，实践中效果也最稳。实现时，Huber鲁棒化实际上等价于给每个量测通道重新分配一个权重——残差小的量测权重接近1，残差大的量测权重被压低。这个权重会进入滤波增益的计算矩阵，从而改变状态更新的方向。
如果你用过加权最小二乘状态估计，会发现这里的思想非常相似。实际上电力系统的静态状态估计里，Huber方法早就被用在了抗差估计中。把它移植到动态状态估计的EKF框架里，逻辑上是通顺的。
2.3 RIEKF的整体滤波流程
RIEKF每个采样周期的计算流程大致分五步：
- 时间更新：用上一时刻的状态和系统方程，预测当前时刻的状态和协方差。
- 量测预测：用预测状态计算量测预测值，并计算雅可比矩阵。
- 残差计算与鲁棒化加权：计算残差，通过Huber函数确定每个量测的权重，构造鲁棒化后的量测噪声协方差矩阵。
- 迭代量测更新：用加权后的量测信息更新状态与协方差，更新后重新线性化，重复多次直到收敛或达到最大迭代次数。
- 输出：取最后一次迭代的状态作为当前时刻的估计值。
从公式层面看，鲁棒化修改的量测噪声矩阵是一个对角矩阵，每个对角元素等于Huber权重与原始噪声方差的乘积。权重计算依赖残差，而残差又依赖当前状态估计，所以迭代和鲁棒化天然是耦合的——这也正是RIEKF把两者结合的深层原因。
3. 电力系统动态状态估计的建模与仿真设计
3.1 系统模型选择：单机无穷大母线还是多机系统
项目里我采用的是经典的单机无穷大（OMIB）系统作为测试平台，发电机用三阶模型描述。三阶模型的状态变量是功角、角速度偏差和暂态电动势，对应的微分方程如下：
- 功角动态方程，反映转子运动；
- 角速度动态方程，包含机械功率与电磁功率的差额；
- 暂态电动势动态方程，反映励磁绕组磁链的变化。
量测配置上，选取发电机功角、角速度和端电压幅值作为量测变量。这套配置跟实际PMU能提供的量测高度一致——功角和频率是PMU最基础的两个输出，电压幅值也是标准量测。仿真时，量测噪声设定为高斯白噪声，同时随机注入若干异常值点，用来检验滤波器的抗差性能。
3.2 场景设计：为什么要加入量测异常值
做对比实验时，很多同学容易犯一个毛病——只在理想噪声环境下比较EKF和RIEKF，结果发现两个滤波器精度差不多，于是得不出主富结论。这个项目我特意设计了三个仿真场景：
- 仅含高斯噪声的正常场景；
- 含个别离群异常值的量测场景；
- 连续多个量测通道同时出现异常值的极端场景。
第一个场景用来验证迭代机制对精度的提升，第二个场景用来验证Huber鲁棒化的抗差效果，第三个场景则检验滤波器在恶劣环境下的稳定性。三组对比下来，EKF、IEKF、RIEKF三者的性能差异一目了然。
从实际工程角度看，PMU数据里的异常值发生频率并没有想象中低。通信延迟、网络拥塞、GPS授时失步，都会导致数据跳变。滤波器不抗差，在线系统就得频繁报警或人工干预，运维成本很高。所以第二个和第三个仿真场景，其实是贴近现场真实情况的。
3.3 评价指标：RMSE与最大绝对误差
对比实验的评价指标我选了两个：全程均方根误差（RMSE）和最大绝对误差（MAE）。RMSE衡量的是整体估计精度，MAE反映的是最坏情况下的偏差——特别是异常值出现瞬间的峰值误差。如果某滤波器RMSE低但MAE很高，说明它在正常情况下表现不错，但遇到异常值会剧烈抖动，这对在线监视系统来说仍是不可接受的。RIEKF的设计目标就是同时压低RMSE和MAE，从实验数据看，效果确实符合预期。
4. Matlab代码实现：从滤波主循环到鲁棒化模块
4.1 代码整体结构规划
写Matlab滤波代码最容易犯的错误是“一把梭”，所有逻辑堆在一个脚本里。项目里我按模块拆分，每个文件只干一件事，调试起来非常省心。代码文件组织如下：
system_model.m：定义系统方程和量测方程；ekf_filter.m：标准EKF实现；iekf_filter.m：迭代EKF实现；riekf_filter.m：鲁棒迭代EKF实现；huber_weight.m：Huber权重计算函数；run_simulation.m：主脚本，负责加载参数、跑仿真、绘图。
如果是在平时的课程作业里，一个脚本从头写到尾倒也跑得通，但只要你后面想换模型、改参数、加对比算法，耦合的代码会让你改一处崩三处。模块化虽然前期多花一点时间，但换来的是整个实验周期的顺畅。
4.2 滤波器主循环的Matlab实现
动态状态估计的主循环本质上是“时间更新—量测更新—数据存储—下一时刻”的循环。RIEKF的主循环里，迭代部分被嵌在量测更新内。下面给出RIEKF核心代码的关键片段：
% 每个采样时刻的滤波循环 for k = 2:T % 时间更新：状态预测与协方差预测 x_pred = f_func(x_est(:, k-1), u(:, k)); A = compute_A(x_est(:, k-1), u(:, k)); % 状态雅可比矩阵 P_pred = A * P_est * A' + Q; % 迭代量测更新的初值设置 x_iter = x_pred; P_iter = P_pred; % 迭代循环：通常3~5次即可收敛 for iter = 1:max_iter % 量测预测与雅可比矩阵 z_pred = h_func(x_iter); H = compute_H(x_iter); % 计算残差并做Huber加权 z_res = z(:, k) - z_pred; weights = huber_weight(z_res, sqrt(diag(R)), 1.345); % 构造鲁棒化的量测噪声矩阵 R_robust = diag(R) ./ weights; % 卡尔曼增益与状态更新 S = H * P_iter * H' + diag(R_robust); K = P_iter * H' / S; x_iter = x_iter + K * z_res; P_iter = (eye(n) - K * H) * P_iter; end % 保存本时刻估计结果 x_est(:, k) = x_iter; P_est = P_iter; end这段代码里的关键操作都集中在迭代循环内部。需要注意的是第10到12行——每次迭代用了新的预测量测和雅可比矩阵，也就是说，迭代一次的“新状态”立刻被用来计算下一迭代的线性化参考点。迭代次数不宜设太大，3到5次通常就足够，过多只会增加计算量，精度提升边际收益很小，甚至可能引入数值振荡。
4.3 Huber权重计算函数的具体实现
Huber权重函数是整个鲁棒化机制的核心，单独拎出来说。函数输入是残差向量、残差标准差和阈值系数，输出是跟残差等长的权重向量：
function w = huber_weight(res, sigma, c) % 标准化残差 r = abs(res ./ sigma); % 分段计算Huber权重 w = zeros(size(r)); idx_small = r <= c; idx_large = r > c; w(idx_small) = 1; w(idx_large) = c ./ r(idx_large); end权重小于1意味着这个量测通道的残差被判定为“异常”，在更新方程里就自动降权。c取1.345是Huber理论的标准建议值，实测里我在PMU数据上试过1.2到1.8之间的不同取值，整体表现差异不大，但1.345确实是兼顾效率和抗差性的甜点区。
很多初学的人会问：为什么不直接把异常量测剔除（置零权重）？我的经验是不要做硬剔除。在动态系统里，一个量测即使残差偏大，也可能携带有效信息，直接扔掉会损失可观测性，尤其是在量测冗余度不高的场景下。Huber降权的好处就是“温柔地忽略”，既能压住异常值的影响，又保留了信息通道。
4.4 参数选取与协方差调参经验
滤波器参数里，过程和量测噪声协方差矩阵的设定对结果影响非常大，而且没有绝对正确的值，只能依据仿真测试反复调试。项目里的调参经验如下：
- 过程噪声Q：状态方程里三个变量的过程噪声设为对角阵，功角的Q值取较小量级，因为功角动态相对平滑；角速度偏差和暂态电动势的Q值可以放宽一点，因为它们受扰动影响更明显。
- 量测噪声R：根据量测通道的信噪比设定，在仿真里用了与真实噪声方差接近的值，以保证滤波器处于信息匹配的状态。R设得过大，滤波器会“不信任”量测，收敛变慢；R设得过小，滤波器会过度信任量测，异常值来临时表现更差。
- 初始协方差P0：给一个适中的对角阵，让滤波器在前几步快速收敛即可。P0设太大会导致最初几个采样时刻的剧烈震荡，太小则收敛慢。
调参的通用原则是：先在无异常值的正常场景下把EKF调解好，再叠加异常值场景调整RIEKF的Huber参数。不建议一上来就直接调RIEKF，否则分不清是哪部分参数导致的性能问题。
5. 仿真结果对比与讨论
5.1 正常场景下：迭代机制带来的精度增益
在仅含高斯噪声的正常场景里，EKF和RIEKF都能稳定跟踪真实状态。两者的RMSE差异并不大，但仔细看每一时刻的估计误差曲线，RIEKF的误差波动明显更小。原因在于迭代更新缓解了强非线性系统里单次线性化的近似误差。
功角估计的RMSE在EKF中大约是0.003 rad量级，RIEKF把这一数值压低了约20%到30%。角速度和暂态电动势的改善幅度略小，这是因为这两个状态的动态方程线性度相对更高，迭代补偿的效果有限。这个结果也印证了一个观点：迭代机制的价值跟系统非线性强度正相关，你的模型越“弯”，迭代带来的收益越明显。
5.2 异常值场景下的抗差能力对比
加入离群异常值后，三类滤波器的表现差异非常突出。EKF在异常值出现的瞬间，功角估计曲线会突然被拉向错误方向，然后花费数个采样周期缓慢恢复；IEKF虽有一定改善，但当异常值残差进入迭代循环后，会被当成有效信息反复利用，收敛速度甚至可能比EKF更慢。
RIEKF的表现则稳定很多。特意在功角量测的第30个采样点注入了一个幅值为真实值两倍的异常值，EKF的估计误差瞬间飙升到接近0.3 rad，而RIEKF几乎看不出明显波动，估计偏差控制在原有噪声水平附近。这就是Huber降权起到的关键作用——异常残差在进入叠代更新前就被压低权重，后续循环里它不会反复捣乱。
如果换成多个量测通道同时出问题，EKF基本会直接发散，而RIEKF仍能维持对状态趋势的跟踪。需要说明的是，RIEKF在极端场景下也会出现误差增大，但增大的幅度和恢复速度远优于其他两类滤波器。从工程角度说，这种“性能软着陆”而不是“直接崩溃”的特性非常有价值。
5.3 鲁棒化与迭代机制的相互作用分析
这套仿真里一个有意思的现象是：迭代机制放大了鲁棒化机制的效果。细想其实不难理解：Huber降权是基于当前迭代点的残差计算的，第一次迭代的权重可能在异常值的“拉扯”下不够准确，但随着迭代推进，状态逐渐回复到合理区域后，残差重新计算、权重随之调整，异常值的影响被逐步剔除。这个协同效应是单纯的鲁棒EKF（不迭代）做不到的。
所以设计RIEKF时，不要试图割裂看待两个机制。你可以做一个消融实验验证这一点：只加Huber不加迭代，异常值来袭时误差峰值会明显压住，但残差尾部会拖得较长；只加迭代不加Huber，正常场景精度提升，但异常值时照样被带偏。两者结合后，既压峰值又缩尾部，才是最终想要的性能曲线。
6. 实操中常见问题与避坑清单
6.1 滤波器发散问题及定位思路
跑仿真遇到滤波器发散，先别急着调参数，按顺序排查下面的可能原因：
第一，检查雅可比矩阵计算是否正确。动态状态估计的雅可比矩阵通常是解析求导推导出来的，手推公式容易出错，尤其是多变量耦合项。可以在代码里加一步数值差分对比，验证解析雅可比和数值雅可比的一致性。这一步能排查掉绝大多数“莫名其妙就发散”的情况。
第二，检查协方差矩阵是否失去正定性。迭代更新过程中，P = (I - K*H)*P这种更新方式在数值上可能产生非对称或非半正定的问题，特别是当量测噪声设置过小、增益过大的时候。建议每步更新后强制对称化处理：P = (P + P')/2，必要时直接改用Joseph形式的协方差更新，数值稳定性会好很多。
第三，检查参数与模型是否匹配。比如仿真用的真实噪声方差是0.01，而你的R矩阵里写的是1，那么滤波器会认为量测非常不可靠，状态更新幅度一直很小，收敛很慢甚至不收敛。反过来也一样，量测噪声设太小，异常值就能肆无忌惮地干扰状态估计。
6.2 迭代不收敛的处理技巧
IEKF和RIEKF在迭代中偶尔会遇到不收敛的问题——迭代若干步后状态估计在某个区间来回震荡，不逐渐逼近稳定点。这通常有两种原因：
一是系统的强非线性导致局部极值点较多，迭代从不同初始点出发可能落入不同吸引域。改善方法是适当降低迭代步长，或者在迭代更新时引入阻尼因子，类似高斯-牛顿法里的阻尼处理。
二是异常值的权重在迭代中切换过猛。比如第1次迭代里某残差恰好落在Huber阈值边缘，权重在1和0.7之间跳变，导致状态在连续迭代中来回抖动。处理方法是在叠代次数不多（比如只剩最后1次）时，将权重固定为第一次迭代计算的结果，避免后期抖振。
6.3 计算效率优化建议
RIEKF的计算量大约是EKF的2到4倍，这在离线仿真里无所谓，但如果要做实时在线应用，需要考虑两个优化点：
- 迭代次数上限设为3即可，第四五次迭代的精度增益通常低于1%，不值得那部分计算开销；
- 雅可比矩阵的复用在迭代后期是非常有效的优化——当相邻两次迭代的状态变化小于预设阈值时，直接沿用上次雅可比矩阵，跳过重算流程。
我在实机上测过，第二个优化几乎无损精度，但能把迭代计算时间压缩掉约30%。如果你的项目对实时性有硬要求，这一步不能省。
7. 扩展到其他场景的思路
RIEKF的适用范围不局限于电力系统。只要你的动态系统满足两个条件——状态方程或量测方程非线性程度较高、量测数据存在被污染的潜在风险——这套“迭代+鲁棒”的框架就能迁移过去。
我实际试过的扩展方向包括：无人机姿态估计中的磁力计异常干扰抑制，把Huber权重的作用对象换成磁场量测残差，效果非常明显；目标跟踪里当目标机动时雷达量测会出现较大残差，RIEKF比标准交互多模型方法更轻量也非常稳定。
在电力系统内部，从单机系统扩展到多机系统并不复杂，只需要把每个发电机的三阶状态变量拼接成一个大状态向量，并在雅可比矩阵里补上互联项。量测配置上如果加入线路潮流量测，滤波器对系统全局状态的感知能力会更完整，代价是矩阵维度和计算量同步上升。从我的实践看，RIEKF在状态维数500以内都能保持良好性能，再往上就需要考虑分布式实现或降维策略了。
如果你手上正好有含异常值标记的PMU实测数据，强烈建议直接把RIEKF拿去测一遍——脱开仿真环境后的表现，会比纯理想场景更有说服力。滤波器实现的本身并不复杂，真正值钱的是面对真实数据时的那点鲁棒性，这也是我做这个项目最深的一个体会。
