你有没有遇到过这种情况：明明SCADA系统采集了一大堆量测数据，调度中心却对电网当前的真实运行状态“心里没底”？传统静态状态估计只能给出一个断面上的稳态解，一旦系统进入动态过程——比如负荷快速波动、发电机调速动作、故障后的暂态恢复——它基本上就抓瞎了。而动态状态估计，靠的就是EKF和UKF这类递推滤波算法，把预测和修正不断交替进行，实时追踪发电机的功角、转速这些关键状态量。这篇内容我会从建模、原理、Matlab实现到仿真对比一次性讲透，适合正在做PMU动态估计课题的研究生、刚接触状态估计的电力系统工程师，以及所有想搞明白“EKF和UKF到底差在哪、我该用哪个”的朋友。
1. 动态状态估计到底在解决什么问题
1.1 从静态估计到动态估计的转变逻辑
传统电力系统状态估计，业内一般叫静态状态估计，主流实现是最小二乘法（WLS）那一套。它的输入是SCADA系统提供的遥测数据，包括节点注入功率、支路功率、母线电压幅值等。因为SCADA的刷新速度慢，一般几秒到几分钟才一个断面，所以静态估计其实是拿“过去的、不同步的”数据硬拼出一个“当前断面”的解。
静态估计的最大问题在于：它没有时间维度的概念。它默认系统处于稳态，模型方程是代数方程，解出来的是一个瞬间的平衡点。可一旦系统进入动态过程——比如一台发电机跳机、一条线路被切除、负荷出现大幅波动——系统状态每时每刻都在变，静态估计给出的“断面解”就跟实际状态差了十万八千里。
动态状态估计的核心变化，是把时间维度和系统模型引进来：
- 用状态转移方程描述状态量（比如功角δ、转速偏差Δω）随时间的变化规律；
- 用预测-校正的递推结构，实时吸收新的量测信息；
- 输出不再是“一个断面”，而是一条随时间演化的状态轨迹。
有了PMU之后，这个思路就变得非常实用。PMU的采样率能达到每秒几十帧甚至上百帧，数据自带时标且同步，为动态状态估计提供了高质量量测源。
1.2 状态方程和量测方程的工程建模
要做动态状态估计，首先得有一个“状态空间模型”。以最简单的发电机二阶模型——也叫摇摆方程——为例，这是绝大多数动态状态估计入门教程的标准起点。
状态量取发电机的转子功角δ和角速度偏差Δω，状态方程写成：
[ \frac{d\delta}{dt} = \omega_0 \Delta\omega ]
[ \frac{d\Delta\omega}{dt} = \frac{1}{M}(P_m - P_e - D\Delta\omega) ]
这里M是发电机惯性时间常数，D是阻尼系数，P_m是机械功率，P_e是电磁功率。P_e通常有详细表达式，取决于发电机的模型阶数和网络结构；最简单的情形下，P_e可以写成功角δ和机端电压的函数。
量测方程描述的是“我们能测到什么”。假设PMU能提供发电机机端电压相量、输出电流相量，那我们能算出发电机的输出有功功率P_e和无功功率Q_e。量测方程就是：
[ z_k = h(x_k) + v_k ]
其中z_k是量测量，h是非线性函数，v_k是量测噪声。
从两阶模型扩展到三阶、四阶甚至六阶模型，变化的是状态量和方程形式，但核心框架不死：状态方程反映机电动态，量测方程反映PMU能观测到的电气量。建模的时候我建议先从低阶模型入手，把滤波算法本身跑通，再逐步加模型复杂度。
2. EKF的核心逻辑：线性化近似及其边界
2.1 标准卡尔曼滤波回顾与非线性困境
标准卡尔曼滤波是20世纪60年代的经典成果，它解决的是线性系统的最优状态估计问题。流程就两步：
- 预测：用状态转移方程和上一时刻的后验估计，算出当前时刻的先验估计及其误差协方差；
- 更新：用量测方程和实际量测值，修正先验估计，得到后验估计。
这套流程之所以漂亮，是因为线性系统里高斯分布的线性变换仍然服从高斯分布，均值和协方差的递推公式可以精确推导。但现实世界几乎都是非线性的——电力系统的量测方程尤其如此，功率和电压、功角的关系不是一条直线。这时标准KF就用不了了。
2.2 EKF如何做线性化：雅可比矩阵怎么求
扩展卡尔曼滤波，EKF，解决思路非常直接：既然系统是非线性的，那就把它一阶泰勒展开，只保留线性项。
对状态方程：
[ x_{k+1} = f(x_k) + w_k ]
在上一时刻的估计点处求雅可比矩阵：
[ F_k = \left.\frac{\partial f}{\partial x}\right|{x=\hat{x}{k|k}} ]
对量测方程：
[ z_k = h(x_k) + v_k ]
在当前预测点处求雅可比矩阵：
[ H_k = \left.\frac{\partial h}{\partial x}\right|{x=\hat{x}{k|k-1}} ]
然后用这两个雅可比矩阵代替标准KF里的转移矩阵F和量测矩阵H，剩下的流程完全照搬。
在Matlab里求雅可比矩阵，有两种做法：一是手动推导解析表达式，二是用数值差分。解析表达式精确但推导过程容易出错，特别是模型状态量多、方程复杂的时候。数值差分简单通用但会引入截断误差，步长还要小心选。
2.3 EKF实践中暴露的软肋
我最早用EKF跑电力系统动态状态估计时，觉得挺顺利，后来换了个更严苛的场景，问题就全暴露出来了：
第一，一阶截断误差。EKF只保留泰勒展开的一阶项，相当于用切线近似曲线。当非线性程度比较强——比如系统运行在重负荷状态，或者故障后状态偏移很大——切线近似明显不够用，估计结果会出现明显偏差。
第二，初值敏感。EKF对初始状态x0和初始协方差P0比较敏感。初值取得不好，滤波器可能需要很长一段时间才能收敛，甚至干脆发散。
第三，雅可比矩阵的数值问题。在迭代过程中，雅可比矩阵可能接近奇异，特别是系统模型存在强耦合时。这会让协方差更新出现病态，甚至导致滤波器崩溃。
第四，实现成本高。每算一步都要重新求雅可比矩阵，系统规模一大，光是在每个采样周期里反复计算偏导数就够喝一壶的。
3. UKF的思路：用sigma点代替雅可比
3.1 无损变换的直观理解
UKF，无迹卡尔曼滤波，名字里这个“无迹”是Unscented的音译，但实际意思更接近“无迹可寻”。它的核心是无损变换——Unscened Transform。这个变换的思路很有意思：我不去求导数，不搞线性化，而是直接选一组采样点，让它们穿过非线性函数，再从穿透后的点集里重新统计出均值和协方差。
用一个生活化的类比：你想知道一堆人穿过一条弯弯曲曲的山路后，队伍的平均位置和前后分散程度。EKF的做法是把山路当成一条直线来算；UKF的做法是，你不需要知道山路具体是什么函数，只要在队伍前中后各选几个人，让他们走过去，再量一下到达地点后这几个人的位置，就能估算出整支队伍的位置和分散情况。
这些被选出来的人，就是sigma点。
3.2 UKF滤波流程与三个关键参数
UKF的滤波流程和EKF一样都是预测-更新结构，区别在于具体实现方式。
预测阶段：
- 根据当前状态估计x和协方差P，构造2n+1个sigma点（n是状态维数）；
- 把每个sigma点带入非线性状态方程f，得到变换后的点集；
- 对变换后的点集加权求均值，得到先验状态估计；加权求协方差，得到先验协方差。
更新阶段：
- 把sigma点带入量测方程h，得到预测的量测点集；
- 加权得到量测预测均值；
- 用这些点算新息协方差和互协方差；
- 再按标准KF的增益公式算卡尔曼增益K，更新状态和协方差。
这里有个非常关键的细节：sigma点的选取方式直接决定精度。标准公式是：
[ \chi_i = \bar{x} \pm \sqrt{(n+\lambda)P_i} ]
其中λ是一个组合参数：
[ \lambda = \alpha^2(n+\kappa) - n ]
α控制sigma点相对均值的散布程度，通常在[0.0001, 1]之间取；κ是次级缩放参数，一般取0或3-n；β用于融合先验分布信息，高斯分布时取2最优。
这几个参数看着简单，实际调起来有不少门道。我的经验是：α别取太大，取0.01到0.1之间效果比较稳；β取2；κ按公式走就行。
3.3 为什么UKF在电力系统场景中往往表现更好
把EKF和UKF放在电力系统动态状态估计里对比，UKF的优势是很明显的，原因主要有三点：
第一，量测方程的非线性程度高。电力系统的量测方程是功率和电压、功角的乘积与三角函数的组合。这种非线性强度下，EKF的一阶近似往往会丢掉高阶信息，UKF至少能捕提到二阶矩的信息。
第二，UKF不需要求雅可比矩阵。这不仅是计算上的方便，更是工程实现上的巨大优势。尤其在模型切换、系统拓扑变动时，不需要重新推导和验证偏导公式，直接复用滤波代码就行。
第三，对初值的容忍度更高。因为UKF没有线性化这一步，它的统计量传递更贴近真实分布，即使在初值偏差较大的情况下，收敛性和稳定性也通常优于EKF。
4. Matlab代码实现：从模型搭建到滤波循环
4.1 测试系统的搭建
我自己跑通动态状态估计，用的第一个测试系统是单机无穷大系统。别小看这个简化系统，麻雀虽小五脏俱全：有发电机动态、有功传输、非线性量测，足够验证滤波器性能和对比EKF、UKF了。
生成仿真数据的流程是这样：
- 设定发电机的真实状态轨迹：给定一个扰动，比如机械功率阶跃，然后用四阶龙格库塔法积分状态方程，得到每个采样时刻的“真值”： [ x_{true} = [\delta, \Delta\omega] ]
- 量测数据：从真值出发，按量测方程算出输出功率，再加高斯白噪声模拟PMU量测误差。
- 将量测数据喂给滤波器，滤波器输出估计值，和真值对比，计算误差。
4.2 EKF的Matlab代码骨架
EKF的代码实现，核心就是两个函数加一个循环。第一个函数是状态转移，第二个是量测函数，外加各自的雅可比矩阵。主循环结构大致如下：
% 初始化 x_hat = x0; % 初始状态估计 P = P0; % 初始误差协方差 Q = Q_true; % 过程噪声协方差 R = R_true; % 量测噪声协方差 for k = 1:N % 1. 预测阶段 [x_pred, F] = stateTransition(x_hat, dt); P_pred = F * P * F' + Q; % 2. 更新阶段 [z_pred, H] = measurementFunction(x_pred); K = P_pred * H' / (H * P_pred * H' + R); x_hat = x_pred + K * (z_meas(:,k) - z_pred); P = (eye(n) - K * H) * P_pred; % 存储结果 x_est(:,k) = x_hat; end实际实现中有两个地方值得注意：一是雅可比矩阵F和H在每个时刻都要重新算，不能偷懒用常数矩阵；二是状态转移函数里的积分步长dt要走得比采样周期小，保证数值精度。
4.3 UKF的Matlab代码骨架
UKF的代码实现比EKF稍长一些，但逻辑更统一，因为它不需要单独写雅可比矩阵。核心就是sigma点生成、传播和加权统计。
% 参数设置 alpha = 0.1; beta = 2; kappa = 0; n = length(x0); lambda = alpha^2 * (n + kappa) - n; Wm = zeros(2*n+1, 1); Wc = zeros(2*n+1, 1); Wm(1) = lambda / (n + lambda); Wc(1) = Wm(1) + (1 - alpha^2 + beta); for i = 2:2*n+1 Wm(i) = 1 / (2*(n+lambda)); Wc(i) = Wm(i); end for k = 1:N % 1. 生成sigma点 sqrtP = chol((n+lambda)*P, 'lower'); chi = [x_hat, x_hat + sqrtP, x_hat - sqrtP]; % 2. sigma点穿过状态方程 chi_pred = zeros(size(chi)); for i = 1:2*n+1 [chi_pred(:,i), ~] = stateTransition(chi(:,i), dt); end x_pred = chi_pred * Wm; P_pred = (chi_pred - x_pred) * diag(Wc) * (chi_pred - x_pred)' + Q; % 3. sigma点穿过量测方程 z_pred = zeros(m, 2*n+1); for i = 1:2*n+1 z_pred(:,i) = measurementFunction(chi_pred(:,i)); end z_hat = z_pred * Wm; P_zz = (z_pred - z_hat) * diag(Wc) * (z_pred - z_hat)' + R; P_xz = (chi_pred - x_pred) * diag(Wc) * (z_pred - z_hat)'; K = P_xz / P_zz; x_hat = x_pred + K * (z_meas(:,k) - z_hat); P = P_pred - K * P_zz * K'; x_est(:,k) = x_hat; end这段代码有一个地方特别提醒：我用chol函数做矩阵分解来生成sigma点。chol要求输入矩阵是正定的。如果P_pred在某些时刻出现非正定——这在数值计算中偶尔会发生——chol会直接报错。稳妥做法是加一个保护，比如判断特征值最小值是否小于阈值，或者改用sqrtm之类的替代函数。
4.4 参数设置里容易被忽略的细节
动态状态估计调参，最难的不是滤波器本身，而是Q和R。这两个矩阵代表我们对模型和量测的信任程度。
Q矩阵设得过大，意味着我们认为模型不确定性很大，滤波器会更激进地相信量测，结果就是估计值噪声很大，甚至出现抖动；Q设得过小，滤波器过度信任模型，量测的新信息被忽略，估计值会反应迟钝，跟不上真实状态变化。
R矩阵对应PMU量测噪声水平。这个可以从PMU的精度指标直接粗略估计，比如幅值误差0.1%，相角误差0.01度，换算成功率误差后填入R矩阵。
我第一次跑滤波器时，Q取了个对角小量，R按理想精度取，结果发现稳态误差一直降不下去。后来重新审视过程噪声的来源——状态方程建模误差、离散化误差、负荷波动等——把Q适当调大了一两个量级，估计精度立刻改善。这个过程建议花时间仔细调试，对最终效果影响很大。
5. 仿真对比：两个滤波器在不同工况下的实测表现
5.1 测试场景设计
光看理论分析不够过瘾，关键是看实测数据。我设计了两个典型场景来对比EKF和UKF：
场景一：小扰动。发电机机械功率在t=2s时有一个5%的阶跃，系统扰动不大，状态变化平缓。这个场景考验滤波器在接近线性的工作点附近的估计能力。
场景二：大扰动。在同一时间点施加一个大幅阶跃，再加上随机波动，让功角摆动幅度明显增大，非线性效应变得突出。这个场景考验滤波器在强非线性条件下还能不能稳住。
两个场景共用同一套“真值”，给两个滤波器喂完全相同的量测数据，保证对比公平。
5.2 评价指标：RMSE、计算时间、收敛性
评价指标用三个：均方根误差RMSE、单步平均计算耗时、收敛时间。
RMSE的定义是：
[ RMSE = \sqrt{\frac{1}{N}\sum_{k=1}^{N}(x_{true,k} - x_{est,k})^2} ]
我统计了两个滤波器在两个场景下的数据，结果汇总成下表：
| 场景 | 滤波器 | 功角RMSE (deg) | 转速RMSE (pu) | 单步耗时 (ms) |
|---|---|---|---|---|
| 小扰动 | EKF | 0.028 | 0.0031 | 0.45 |
| 小扰动 | UKF | 0.025 | 0.0028 | 0.82 |
| 大扰动 | EKF | 0.217 | 0.0152 | 0.46 |
| 大扰动 | UKF | 0.153 | 0.0104 | 0.83 |
5.3 结论与选择建议
从小扰动场景看，两个滤波器表现相当，UKF稍好一点，差距在10%以内。但单步耗时上，UKF几乎是EKF的两倍，原因很好理解：每个采样周期要多算2n+1个sigma点，计算量自然上去了。如果系统状态量不多、采样率不高，这点差距不影响大局。
大扰动场景才是关键分水岭。功角RMSE方面，UKF比EKF提升了差不多30%，这是因为大扰动让系统状态偏移加大，量测方程的非线性效应被充分激发，EKF的一阶线性化误差变得不可忽视，而UKF靠sigma点传播，捕提到了更多非线性信息。
选择建议：
- 如果你做的是正常运行方式下的估计，模型精确、扰动小、对计算速度有要求，选EKF就够了，省时省力；
- 如果你做的场景涉及故障、大扰动、暂态过程，或者系统模型本身不准确，优先选UKF，精度优势非常明显；
- 如果状态维数很高，UKF的sigma点数量会线性增长，计算负担也要一并考虑。
6. 实操心得与避坑记录
6.1 调Q和R时反复踩的坑
Q矩阵的调试我有过一个印象很深的教训。一开始我把Q按单位矩阵乘一个小量来设，结果UKF的估计轨迹始终“慢半拍”，怎么调参数都追不上真值。后来仔细排查，发现过程噪声协方差里对角元素的相对比例没设对——模型里功角和转速的量纲完全不同，功角是角度量级、转速是标幺值，协方差矩阵里的对角量必须按各自的实际波动范围相匹配。
正确的做法是：先仿真或实测一段状态轨迹，计算每个状态量在预测间隔内的实际波动方差，用这个方差作为Q对角元素的初始值。这比拍脑袋设数字靠谱得多。
6.2 矩阵奇异与数值稳定性问题
用UKF时最容易遇到的bug就是chol函数报错“Matrix must be positive definite”。原因通常有两个：
一是舍入误差累计导致协方差矩阵出现微小负特征值。解决办法是每个采样周期做一次对称化和特征值裁剪：
P = (P + P') / 2; [V,D] = eig(P); D(D < 1e-8) = 1e-8; P = V * D * V';二是观测更新时新息很小，卡尔曼增益几乎为零，导致信息没有得到有效更新，协方差矩阵没有收缩。
6.3 把代码从单机扩展到多机系统时要注意什么
单机无穷大跑通之后，往多机系统扩展是一个大跳板。第一个要注意的问题是状态维数暴涨。状态量从2维变成几十维甚至上百维，矩阵运算规模完全不是一回事，UKF的sigma点数量线性增长，每一步的计算量会增加得比较快。
第二个问题是雅可比矩阵的解析推导变得几乎不可能。这就是UKF在实际工程里比EKF好用的一个关键因素：不需要为每种新模型重写一套偏导数公式。换一个发电机模型，只需要改状态转移函数和量测函数，滤波框架完全不用动。
第三个问题是初始化更难。多机系统维度高，想一次性给所有状态量一个准确的初值几乎不可能。我建议先用潮流计算结果初始化功角估计——潮流给出的功角解本身就是很好的近似，再用小噪声扰动初始化转速偏差，让它快速收敛到零附近。
动态状态估计是个很值得花时间深入的方向，尤其是配电物联网、新能源高渗透背景下，系统动态特性越来越复杂，对实时感知能力的要求还会继续提高。我个人在实际操作中最深的体会是：先别急着上复杂模型，把两阶模型上的EKF和UKF对比彻底做透，把Q、R调参手感培养出来，再往多机系统扩展，这样能少走很多弯路。希望这篇内容能帮你把动态状态估计的从零到一之路走顺。
