做电力系统动态状态估计，EKF和UKF是两个绕不开的经典选项。我最近正好用Matlab完整实现了一套基于扩展卡尔曼滤波（EKF）和无迹卡尔曼滤波（UKF）的电力系统动态状态估计器，把单机无穷大系统的机电暂态过程作为测试案例，从模型搭建、量测生成到滤波器闭环全部走了一遍。这篇博客就把这套实现拆开讲清楚：动态状态估计到底在估什么、EKF和UKF的原理差异、Matlab代码结构和参数调试方法，以及我实际调试过程中踩过的坑。适合正在做PMU量测数据处理、机电暂态辨识，或者拿动态状态估计做课程设计的研究生和工程师参考。
1. 项目到底在做什么：动态状态估计的定位与核心难点
1.1 什么是动态状态估计，和静态估计差在哪
传统电力系统状态估计大多基于加权最小二乘（WLS），它假设系统处于稳态，只对某一个时间断面的SCADA量测做处理，输出的是一组静态潮流解。静态估计本身没有时间维度，状态量之间没有动态关联，量测少了或者坏数据多了，只会影响当前断面的精度。
动态状态估计则完全不同。它把发电机转子运动、励磁调节等动态过程写进状态方程，状态量通常是功角、角速度、暂态电动势、励磁电压等，然后用量测数据不断修正预测结果。说直白一点，静态估计像拍照，动态估计像拍视频。二者最大的差异在于动态估计引入了“状态随时间演化的规律”，也就是系统模型本身，这让它在暂态过程中依然能给出平滑、可信的轨迹。
我这次项目的测试场景是单机无穷大系统，动态模型采用三阶实用模型。状态量取发电机的功角、角速度偏差和q轴暂态电动势，量测来自模拟的PMU数据。之所以选这个规模，是因为单机场景所有物理量都容易推导，滤波器好坏一眼就能看出来，调试阶段不会因为模型庞大而分不清是建模问题还是滤波问题。
1.2 为什么同时实现EKF和UKF，而不是只做一个
很多教材只讲EKF，因为EKF是卡尔曼滤波在非线性系统上的最直接扩展，思路就是把非线性的状态方程和量测方程在估计点附近做一阶泰勒展开，然后用线性卡尔曼滤波的框架。它的代码量小，计算量低，控制工程里用了很多年。
但电力系统发电机模型尤其是功角、有功功率之间的关系，非线性程度相当强。EKF用切线代替曲面，模型本身误差就摆在那里，一旦系统进入深饱和或者大扰动场景，线性化误差会被协方差矩阵的传播进一步放大，滤波结果可能出现明显偏差甚至发散。
UKF走的是另一条路：不计算雅可比矩阵，而是用一组精心构造的Sigma点去直接传播均值和协方差。它的思想是用样本点逼近非线性函数的统计特性，而不是用一阶导数逼近函数本身。这种方法在三阶截断精度上通常优于EKF，尤其在强非线性场景下优势更明显，代价是计算量更大。
所以这个项目同时实现两个滤波器的意义在于：在同一个模型、同一组量测、同一套噪声矩阵下做公平对比，既能验证UKF在精度上的优势，也能量化它比EKF多花的计算时间。对于实际工程来说，这种对比很重要，因为在线动态状态估计往往有时间窗限制，不是精度越高越好，而是要在硬件资源允许的范围内选最合适的滤波器。
2. 数学模型与滤波原理：公式背后的物理意义
2.1 发电机组动态模型怎么建
动态状态估计第一步是把物理系统写成离散状态空间形式。这里以三阶单机无穷大模型为例，标幺值下的状态向量为：
x = [δ, ω, Eq′]^T
其中δ是发电机功角，ω是角速度（标幺值，基准为同步转速），Eq′是q轴暂态电动势。连续时间状态方程可以写成：
dδ/dt = ω0 * (ω - 1)
dω/dt = (Pm - Pe - D * (ω - 1)) / (2H)
dEq′/dt = (Efd - Eq′) / Td0′
这里ω0是同步角速度，Pm是机械功率，Pe是电磁功率，D是阻尼系数，H是惯性时间常数，Efd是励磁电压，Td0′是励磁绕组暂态时间常数。电磁功率Pe的表达式为：
Pe = Eq′ * V / Xd′ * sin(δ)
其中V是无穷大母线电压，Xd′是暂态电抗。这是经典的三阶模型，做学术研究和初阶工程验证完全够用。
量测方程方面，PMU能够给出的典型量测包括功角、电压幅值、有功功率和无功功率。为了便于观察，我这里的量测向量取为：
z = [δ, Vt, Pe]^T
其中Vt是机端电压幅值，它和Eq′、功角δ之间存在代数关系。把状态方程和量测方程写成离散形式时，需要选择一个采样周期Δt，用前向欧拉或四阶Runge-Kutta做离散化。我实测下来，前向欧拉在采样周期小于1毫秒时精度足够，但为了贴近真实PMU数据（常见上报率有10ms、20ms），建议用四阶Runge-Kutta生成真值轨迹，同时滤波器内部预测时也采用同样的离散化方式，避免数值积分误差被滤波器误判成过程噪声。
2.2 EKF核心逻辑：线性化像用切线近似曲面
EKF的完整流程分为预测和更新两步。
预测步用上一时刻的估计值通过状态方程得到先验估计：
x̂_k|k-1 = f(x̂_k-1)
P_k|k-1 = F_k-1 * P_k-1 * F_k-1^T + Q
其中F_k-1是状态方程对状态向量求偏导得到的雅可比矩阵，Q是过程噪声协方差矩阵。这一步的本质是：用模型“猜”下一时刻状态，同时把估计的不确定性也往前推一步。
更新步引入量测：
K_k = P_k|k-1 * H_k^T * (H_k * P_k|k-1 * H_k^T + R)^(-1)
x̂_k = x̂_k|k-1 + K_k * (z_k - h(x̂_k|k-1))
P_k = (I - K_k * H_k) * P_k|k-1
其中H_k是量测方程对状态的雅可比矩阵，R是量测噪声协方差矩阵。卡尔曼增益K决定了滤波器相信模型还是相信量测：R相对Q越大，K越小，滤波结果越平滑；反之K越大，滤波结果越容易跟随量测噪声。
在Matlab实现EKF时，最麻烦的是F和H的解析推导。我的做法是用符号工具箱Symbolic Math Toolbox生成雅可比矩阵，然后转换成普通的matlabFunction，这样既能避免手推出错，又能在参数变化时快速重新生成。如果不用符号工具箱，也可以采用数值差分法：对每个状态分量加一个小扰动，重新计算状态方程和量测方程的差值。差分法鲁棒性尚可，但扰动步长要反复试，太大会引入截断误差，太小会引入舍入误差，建议优先用符号求导。
2.3 UKF核心逻辑：用Sigma点群代替线性化
UKF的思路是：既然状态分布未知，干脆用一组确定的采样点去“代表”这个分布，让这些点分别通过非线性函数，再统计输出的均值和协方差。对于n维状态向量，需要2n+1个Sigma点。
Sigma点选取的经典公式为：
χ0 = x̂
χi = x̂ + (sqrt((n + λ) * P))_i, i = 1,...,n
χ(n+i) = x̂ - (sqrt((n + λ) * P))_i, i = 1,...,n
这里的(sqrt(...))_i表示矩阵平方根的第i列。λ是复合缩放参数，通常定义为λ = α^2 * (n + κ) - n。α是一个很小的正数，控制Sigma点离均值有多远，一般取1e-3到1量级；κ一般取0或3-n；β用来融合先验分布信息，高斯分布时取2。
每个Sigma点都有一对均值权重Wm和协方差权重Wc，具体权重公式可以在任意一本讲UKF的书里查到。在Matlab中计算矩阵平方根我推荐用chol函数做Cholesky分解，或者用sqrtm计算通用矩阵平方根。有一点必须小心：P矩阵必须是对称正定的，否则算出来的Sigma点会明显偏移，导致滤波发散。
所有Sigma点生成后，将每个点代入状态方程得到一组预测点，然后加权求均值和协方差。更新步同样将每个预测点代入量测方程求取量测预测点，再计算协方差和互协方差，最终求出卡尔曼增益。这个过程中看不到雅可比矩阵，但它的数学本质是“用采样近似分布”，当Sigma点数量和分布参数选得恰当，得到的均值和协方差精度可以达到三阶泰勒展开，比EKF的一阶线性化高一个量级。
我在实现UKF时踩过一个具体坑：如果过程噪声Q取得特别小，预测协方差矩阵P会出现特征值接近零甚至为负的情况，Cholesky分解直接报错。后面我会单独说怎么应对。
2.4 EKF与UKF精度和计算量对比
从理论上说，EKF适合非线性较弱的场合，UKF适合强非线性场合。但在电力系统动态状态估计这个场景里，我实测的几个结论值得参考：
| 对比维度 | EKF | UKF |
|---|---|---|
| 实现复杂度 | 低，代码量小 | 中等，需要处理Sigma点 |
| 模型导数依赖 | 需要解析或数值Jacobian | 完全不需要 |
| 三阶以上非线性 | 难以处理 | 自然包含 |
| 单步计算时间 | 快 | 大约慢3-5倍 |
| 协方差保持 | 相对稳定 | 需要正定约束 |
| 强暂态场景 | 可能偏差明显 | 明显更稳 |
| 适合场景 | 实时性要求高、模型弱非线性 | 精度要求高、算力充足 |
这个表格不是绝对真理，但它能帮你做第一步决策。如果目标是水电机组、火电机组的在线动态仿真，采样频率非常高且硬件算力受限，EKF往往更现实。如果是离线分析故障录波，或者做算法对比研究，UKF的精度优势更容易体现出来。
3. Matlab 代码实现细节：从数据生成到滤波闭环
3.1 整体文件结构和接口设计
我习惯把整个项目拆成清晰的模块，而不是把所有代码塞进一个main脚本里。仿真工程结构大致如下：
project/ main.m % 主程序：设置参数，生成量测，运行滤波器，绘图 system_model.m % 状态方程和量测方程（连续/离散） gen_true_trajectory.m % 用ode45生成真值状态轨迹 gen_measurements.m % 由真值加噪声生成量测 ekf_filter.m % EKF滤波器 ukf_filter.m % UKF滤波器 plot_results.m % 绘图和误差分析关键设计点是EKF和UKF的接口必须完全一致。我定义统一调用格式：
[x_est, P_est, stats] = ekf_filter(x0, P0, Q, R, state_fun, meas_fun, z, dt); [x_est, P_est, stats] = ukf_filter(x0, P0, Q, R, state_fun, meas_fun, z, dt);其中state_fun和meas_fun都是函数句柄，分别代表离散状态转移和量测计算。这样切换滤波器只需要改一行调用语句，所有评估和绘图代码都不用动。
两个滤波器内部都维护一个时间序列结构体，记录每一时刻的状态估计、协方差、卡尔曼增益、残差、计算耗时。这个信息后面做收敛性诊断非常有用，不要只存状态估计值。
3.2 初始化、噪声矩阵和参数调整
初始化是动态状态估计最容易做错的环节。P0代表初始状态估计的可信程度，如果P0设置得过小，滤波器会“懒得”纠正初始误差，导致前几十个点明显偏离真值。我一般把P0设成对角阵，对角线元素根据状态量纲取不同的量级。功角单位是弧度，初始不确定度可以给0.01到0.1；角速度标幺值的不确定度给1e-4到1e-3；暂态电动势给0.01左右。
Q矩阵的设置更要小心。Q代表过程噪声，也就是“模型本身有多不准确”。如果系统模型很准，Q取小值，滤波轨迹很平滑；如果模型有模型误差，Q取大值，滤波器会更相信量测。实际项目里Q常常不是简单对角阵，因为状态方程之间相互耦合，比如机械功率突增会导致功角和角速度同时变化，这种耦合可以体现在Q的非对角元素上。但初始调试时建议先设对角阵，把对角线逐项调整到合理范围，再考虑耦合项。
R矩阵相对直观，因为它对应PMU量测误差。功角量测误差协方差可取(0.001 rad)^2，电压幅值取(0.001 p.u.)^2，有功取(0.01 p.u.)^2。如果做仿真时对量测噪声添加了不同强度，R必须和实际添加噪声的方差匹配，否则滤波器的增益计算就是错的。
一个我反复用的小技巧：先把R固定成已知量测方差，把Q从非常小的值开始逐渐增大，画出状态估计的残差曲线。残差如果表现出明显的低频偏移，说明Q太大导致滤波器过度相信量测；如果残差白噪声但状态曲线滞后，说明Q太小，滤波器更新太慢。这个经验在调试任何卡尔曼类滤波器时都适用。
3.3 仿真场景：短路故障与评估指标
单机无穷大系统是最简单的动态估计测试场景，但简单不意味着没有挑战。我设计了一个典型场景：系统在t=0.5s发生三相短路故障，持续0.1s后切除，整个机电暂态过程持续5秒。故障期间电磁功率变化剧烈，功角和角速度快速震荡，这正好考验滤波器在强非线性过程中的跟踪能力。
生成真值轨迹用ode45，量测则按10ms采样率在真值上加高斯白噪声。滤波器同样按10ms步长运行。这里有个细节：量测向量可能不是每一拍都完整，实际PMU偶尔会有丢帧，我特意在部分时间点把量测向量置为NaN，然后在滤波器中做“无量测更新”处理。这个功能很实用，后面排查时也有用。
评估指标我用三个：
- RMSE：状态量估计值与真值之间均方根误差，衡量整体精度。
- 最大绝对误差：观察最坏情况下滤波器的偏离程度。
- 单步平均耗时：用tic/toc统计每步滤波所需时间，对比EKF和UKF的计算速度。
这三个指标不需要特别复杂的代码，但用于横向对比滤波器已经足够。画图时我分成上下两排：上排画真值、EKF、UKF的状态曲线，下排画误差绝对值曲线。实际效果通常是：EKF在故障瞬间误差瞬时冲得比较高，UKF相对平缓；收敛之后两者误差都会回落，但UKF稳态误差更小。
4. 踩坑记录：EKF/UKF实现中容易翻车的地方
4.1 协方差矩阵发散：现象、原因和解决
我在第一次跑UKF时遇到最经典的问题：运行到故障切除瞬间，P矩阵出现负特征值，Cholesky分解报错。后来定位到两个原因：一是故障瞬间状态突变剧烈，模型预测误差远超Q的假设；二是数值积分过程把协方差矩阵的对称性破坏了。
解决办法有三个。第一，每次更新后强制P矩阵对称化，即P = (P + P^T) / 2。第二，在Cholesky分解前加一个小对角扰动，比如P = P + 1e-12 * eye(n)，避免特征值刚好为零。第三，如果P矩阵实在不正定，就采用协方差的平方根形式，这就是平方根UKF（SRUKF）的基本思路。虽然本项目没有用到SRUKF，但理解这一点能帮助排查类似问题。
EKF也会遇到发散问题，表现是残差越来越大，状态曲线直接飞出物理范围。此时优先检查Q是否太小，以及初值P0是否和实际初始误差匹配。如果调试时让滤波器从错误初值启动，P0又给得特别小，滤波器会僵住，纠正不过来，这是初学者最容易踩的坑。
4.2 雅可比矩阵算错导致EKF结果诡异
EKF的雅可比矩阵是手推最容易出错的地方，尤其是量测方程里同时有sin(δ)和Eq′乘积的偏导数。我对EKF的建议是：不要手算，直接用Matlab Symbolic算一遍，并且把符号表达式打印出来和手推结果校验一次。
还有一种常见的工程做法是数值差分求雅可比。核心代码类似于：
F = zeros(n, n); for j = 1:n xp = x; xp(j) = xp(j) + epsJ; xm = x; xm(j) = xm(j) - epsJ; F(:, j) = (state_fun(xp) - state_fun(xm)) / (2 * epsJ); end扰动步长epsJ根据状态量取值范围取相对小量，比如1e-6 * max(abs(x))。这个方法的优点是通用，缺点是计算量翻倍，而且步长选不好会让雅可比矩阵出现噪声，导致增益抖动。我最终还是在EKF里用了符号求导生成的函数，运行速度和稳定性都更好。
如果你和我一样做对比实验，还要注意一个问题：EKF的雅可比矩阵每一步都需要在当前估计点重新计算，不能复用上一时刻的矩阵。故障切除瞬间状态跳变很大，用旧矩阵会导致预测方差完全不准。
4.3 UKF的Sigma点权重和正定性
UKF的两个常见“暗坑”是权重归一化和Sigma点对称性。权重公式本身已经满足归一化，但如果你修改了α、κ参数，最好在初始化时显式计算一次权重总和，输出到命令行检查是否等于1。我见过有人在代码里忘加负的权重项，导致均值和协方差估计全部偏移。
另外，Sigma点必须严格关于均值对称。如果P的平方根计算有数值误差，Sigma点会不对称，滤波结果随之偏差。每次生成Sigma点之后，我习惯用mean(sigma,2)和x̂做比较，如果差异超过1e-8，就说明P矩阵处理有问题。
UKF还有一个容易被忽略的点：过程噪声的注入位置。正确做法是在状态方程传播Sigma点之后再添加Q，也就是预测协方差里加上Q。如果Q重复加了两次，或者忘了加，滤波器一个会过于平滑，一个会过于发散。这些细节不写一遍代码根本发现不了。
4.4 采样步长和量测频率不一致
实际PMU数据上报频率和状态方程的自然时间尺度往往不一致，比如状态方程要用1ms步长做数值积分，但PMU只能每20ms给一个量测点。此时正确的做法是：在两次量测之间执行多次预测步，只有拿到量测时才执行更新步。这个概念叫多步预测单步更新。
我最初实现时图省事，直接把滤波步长设成和量测周期一样，结果强非线性场景下数值积分误差明显，状态曲线出现锯齿。后来改进为：滤波器内部用固定小步长进行预测，量测到达时更新时间，效果立刻改善。
代码实现上也不复杂：在循环里判断当前时刻是否有量测，有则调用更新函数，没有则跳过更新。这个功能建议一开始就写在滤波器接口里，不然后面做真实PMU数据处理时还得重构。
4.5 问题速查表
我把调试过程中遇到的高频问题汇总成了一张表，方便你对照排查：
| 现象 | 可能原因 | 快速处理 |
|---|---|---|
| 滤波曲线严重滞后 | Q太小，或R太大 | 减小R，增大Q |
| 滤波曲线抖动剧烈 | R太小，或状态初值不准 | 增大R，重新设置P0 |
| P矩阵报非正定 | 数值对称性被破坏 | 强制P=(P+P')/2，加对角扰动 |
| EKF误差比UKF大很多 | 模型非线性强 | 检查雅可比是否在最新点计算 |
| UKF耗时超标 | Sigma点数量多 | 尝试减少状态维数，或考虑SRUKF |
| 故障瞬间误差冲高 | 模型突变提前被Q低估 | 适当增大Q或引入模型切换 |
这个表看起来简单，但里面每条都是我实际跑Matlab时遇到并解决的，不是抄文档。
5. 从单机到多机：一点经验与扩展建议
5.1 P0、Q和R调整的实战经验法则
动态状态估计的精度有上限，这个上限由“过程噪声Q、量测噪声R、模型准确度”三者共同决定。经验法则是：先在离线仿真里把R设定为已知量测方差，然后将Q从1e-6量级逐渐增大，画出RMSE关于Q的曲线。你会发现RMSE随Q增大先下降后上升，存在一个明显的“盆地”，最优Q通常在这个盆地底部附近。这个操作听起来简单，但对每个状态维度都要分别做一遍，统一取一个Q矩阵反而容易顾此失彼。
还有一个被低估的参数是P0。如果滤波器从很差的初值启动，P0给太大会导致前几十个点误差很大，给太小又纠正不过来。我的建议是P0按照初值偏差平方的2到3倍来设定。比如初始功角估计可能偏差0.1 rad，那P0(1,1)就设置在0.02到0.03。
5.2 扩展到多机系统与改进滤波
单机模型调试没问题后，往多机系统扩展会遇到两个主要问题：状态维数上升导致计算成本增加，以及量测方程更复杂导致雅可比矩阵推导更困难。EKF在多机场景的雅可比矩阵可以用稀疏矩阵技术加速，UKF则可以尝试降维采样或者平方根形式。
如果Q和R都无法准确获取，可以考虑自适应滤波，也就是在线估计噪声统计特性。常见做法是使用新息序列（innovation sequence）的协方差来修正R，使用状态残差统计来修正Q。这个思路可以在不改变滤波器主结构的前提下显著提高鲁棒性，但对编码能力要求更高。
在实际工程中，除了EKF和UKF，还可以关注集合卡尔曼滤波（EnKF）、粒子滤波和滚动时域估计（MHE）。它们各有优劣，但EKF和UKF依然是性价比最高的起步方案。作为电力系统方向的研究者，先把这两种掌握扎实，后面再切换到任何高级滤波器都会容易得多。
5.3 残差统计检验还能告诉你什么
滤波器是否匹配，不能只看曲线贴不贴真值。一个更系统的办法是统计新息向量（innovation）的均值、方差，并用卡方检验判断它是否符合理论分布。如果新息均值明显偏离零，说明系统模型或量测模型有偏；如果新息方差远大于理论值，说明R被低估；如果远小于理论值，说明R偏大或者Q偏小。这个诊断方法在写论文和做工程报告时都非常有用。
我用Matlab实现统计检验时，只用了几个标准函数：mean、cov、chi2inv。画出新息序列和卡方检验限带后，能从图上直接看出滤波器某个时间段是否存在异常。这个方法比单纯看RMSE更能解释“为什么会变差”。我在项目的README里专门加了这样一行实验记录，后来回看时发现很多问题都是用这个方法定位的，强烈建议你也养成这个习惯。
这个项目做到最后，我的体会是：EKF和UKF的代码实现并不难，难的是对模型物理过程的理解，以及把数值细节处理到位。如果你只在理想仿真里调通参数，那还不够；建议你试着把量测改成随机丢帧、把初值故意设置偏差、把系统参数人为改成错误值，再观察两种滤波器的表现。这样折腾一轮之后，你对Kalman框架的理解会完全不一样。最后再分享一个小技巧：所有涉及P矩阵的运算都用MEX或者向量化写法，Matlab循环在状态维数提高后会明显拖慢UKF的速度，提前优化能让你后续扩展少走很多弯路。
