扩展卡尔曼与无迹卡尔曼滤波：电力系统动态状态估计实战解析

电力系统状态估计从“静态断面”走向“动态过程”，正在成为调度自动化里越来越绕不开的一项技术。尤其是同步相量量测单元（PMU）普及之后，量测数据的时间分辨率从秒级提升到几十毫秒级，如果仍然用传统的加权最小二乘静态估计，很多暂态过程根本来不及被捕捉，更谈不上给在线决策提供支撑。这篇博文不打算按教科书顺序从头推导卡尔曼滤波公式，而是站在Matlab代码实现的角度，把扩展卡尔曼滤波（EKF）和无迹卡尔曼滤波（UKF）在电力系统动态状态估计中的应用细节拆开讲清楚。你会看到完整的滤波主循环骨架、雅可比矩阵和Sigma点传播的具体写法、算例对比的典型结论，以及我在实际调参过程中踩过的坑。内容适合正在做动态状态估计仿真、课设或论文复现的读者，也适合想从静态潮流计算转向动态分析方向的工程师参考。

1. 传统静态估计的瓶颈与动态状态估计的登场

1.1 电力系统状态估计到底在解决什么问题

状态估计在能量管理系统里承担的角色，相当于给调度员一个“看得懂的全景图”。量测装置分散在各个变电站，但量测值总是带噪声、带坏数据，甚至不同厂站之间的采集时刻都不完全一致。状态估计的任务是把这些零散的带噪声量测，结合网络拓扑和线路参数，计算出全网各节点电压幅值和相角——也就是所谓的“状态向量”——让调度、潮流分析、安全评估都能工作在一个自洽的断面上。

传统做法绝大多数是加权最小二乘类方法：把非线性量测方程迭代线性化，然后求解一个优化问题，得到某一时刻的静态断面。这个思路在断面变化缓慢、量测数据每几分钟才刷新一次的传统数据采集与监控系统时代完全够用。上一轮量测算完得到一组状态，下一轮量测来了再从头算一遍，两道工序之间没有任何时间维度的关联。

1.2 为什么静态方法应付不了动态场景

新能源大规模接入、负荷快速波动、电网运行点频繁变化之后，静态估计的两大局限被放大了。

第一，它没有利用历史信息。每个断面都是独立的优化问题，上一时刻的结果对当前时刻完全没有帮助，相当于每次考试都从零开始做。第二，它对动态过程的描述能力很差。低频振荡或者有功突增之后，发电机的功角和转速会在几十秒内连续摆动，这时候用静态估计只拿某一秒的量测去解一个断面，得到的是“这一刻最可能的稳态值”，而不是系统此刻真实的动态轨迹。

动态状态估计的思路完全不同：它把状态看成随时间演化的变量，用系统的动态模型来描述状态转移，再用新的量测对状态预测做修正。预测由模型给出，更新由量测给出，两者按各自的不确定性加权融合。卡尔曼滤波就是这种“预测+更新”结构最经典的实现，而在电力系统这种强非线性场景下，扩展卡尔曼滤波（EKF）和无迹卡尔曼滤波（UKF）是应用最广、实现成本相对可控的两种方案。

2. 把发电机动态写进状态空间：建模与离散化细节

2.1 状态变量、输入与输出的约定方式

动态状态估计的目标通常是同步发电机的状态。最常用的基础模型是二阶经典摇摆方程：

d(delta_i)/dt = omega_i - omega_s

d(omega_i)/dt = (omega_s / (2H_i)) * (Pm_i - Pe_i - D_i * (omega_i - omega_s))

其中delta_i是第i台发电机的功角，单位弧度；omega_i是电角速度，单位弧度每秒；omega_s是同步转速；H_i是惯性时间常数，单位秒；D_i是阻尼系数；Pm_i是机械功率输入，标幺值；Pe_i是电磁功率，标幺值。

状态向量取所有发电机的功角和转速，维度是发电机数量的两倍。输入u通常取机械功率和励磁电压这类对动态过程起驱动作用的量。量测量z来自同步相量量测，最常见的是母线电压幅值和相角，部分场景下也直接测量发电机内电势或支路功率。

这个二阶模型看起来简单，但足够检验EKF和UKF的算法特性。如果要进一步细化，比如加入暂态电抗后的暂态电动势、q轴和d轴电压分量，甚至接入励磁系统、调速器模型，滤波器维度和非线性程度都会显著上升，但整体框架完全不变。我的建议是先拿经典二阶模型把算法流程跑通，再逐步增加模型复杂度，不要一上来就堆一个高维详细模型，那样出了问题根本分不清是算法问题还是模型问题。

2.2 连续微分方程的离散化处理

卡尔曼滤波工作在离散时间序列上，所以必须把连续微分方程转换成离散的状态转移函数。最简单可行的是欧拉法：

delta_i(k+1) = delta_i(k) + (omega_i(k) - omega_s) * dt

omega_i(k+1) = omega_i(k) + (omega_s / (2H_i)) * (Pm_i - Pe_i(k) - D_i * (omega_i(k) - omega_s)) * dt

dt取决于量测采样间隔。PMU典型输出是每秒30帧到60帧，也就是dt约等于0.033秒或0.017秒。对二阶摇摆模型来说，这个步长在数值上足够稳定。但注意，一旦模型里加入励磁动态，其中一些状态的时间常数较小，欧拉法就有点吃紧了，建议换四阶龙格库塔法。

这里有个新手特别容易踩的坑：欧拉法虽然写法简单，但当阻尼项或输入变化很剧烈时，积分误差会逐步累积到状态预测里，导致滤波器预测误差偏大。卡尔曼更新虽然能部分校正回来，但预测这一步持续输错，滤波器会处于被动状态。比较好的做法是把状态转移函数独立封装成一个函数文件，后续升级成龙格库塔时只改这一个文件，主循环不用动。

2.3 非线性量测方程与雅可比矩阵的准备

量测方程h(x)负责把状态向量映射到量测空间。典型实现里，由发电机功角得到内电势相量，结合网络导纳矩阵求母线电压，再由母线电压计算节点电压幅值、相角以及支路功率。

这个映射高度非线性，具体形式取决于量测装置配置。EKF需要求h对x的雅可比矩阵H，这一步是手推公式最容易出错的地方，漏一个偏导或者符号写反，滤波器直接发散。我自己做Matlab代码验证时一定会用中心差分做数值对照，把数值雅可比和解析雅可比放在一起比较，确认每个元素都一致后才继续跑滤波主循环。

2.4 噪声协方差矩阵的初始设定

过程噪声协方差Q和量测噪声协方差R的初始值，直接影响滤波器的“信任偏好”。PMU的电压幅值量测误差通常在0.1%量级，相角误差在0.01弧度量级，把这些物理量级换算成方差可以直接作为R的对角元素。Q的物理含义是状态转移模型没能描述的那部分动态，比如负荷随机波动带来的机械功率变化。取值过大会让滤波器不信任模型，增益变大，状态轨迹受量测噪声影响明显；取值过小则会让滤波器过于信任模型，发生“锁死”在错误状态上的情况。

这里还要注意单位统一。相角如果是弧度，整个代码里就统一用弧度；如果混用了度和弧度，后面调参的时候你会被看起来很奇怪的结果折腾很久。

3. EKF与UKF的分岔点：线性化误差与无损变换

3.1 EKF的雅可比矩阵到底哪里容易出问题

EKF的核心假设是：非线性系统在当前估计点附近可以近似为线性。具体做法是把状态转移函数和量测函数在当前估计点做一阶泰勒展开，用雅可比矩阵替代线性卡尔曼滤波里的状态转移矩阵A和量测矩阵H。

问题在于电力系统的量测方程非线性很强，暂态过程中运行点变动又非常剧烈，一阶近似带来的截断误差不容忽视。误差的实际影响不只是精度下降，更危险的是滤波器一致性被破坏：线性化点离真实状态较远时，雅可比矩阵失效，卡尔曼增益算得不准确，协方差矩阵一步步变小，滤波器以“高度自信”的姿态锁定在错误状态上，工程上说的滤波器发散往往就是这个过程。

3.2 UKF的Sigma点传播与权重设计

UKF换了一个思路：不求解析雅可比，而是用一组精心选取的采样点——Sigma点——来刻画当前状态分布。把这组点分别通过非线性函数传播，再用传播后得到的点集重新拟合均值和协方差。

Sigma点的生成方式如下。设状态维度为n，选择参数alpha、beta、kappa，计算：

lambda = alpha^2 * (n + kappa) - n

然后在状态均值附近取2n+1个采样点。第一个点就是均值本身，其余2n个点在均值两侧，偏移量由矩阵(n+lambda)*P的Cholesky分解列向量确定。

权重分配是：

Wm_0 = lambda / (n + lambda)

Wc_0 = Wm_0 + (1 - alpha^2 + beta)

Wm_i = Wc_i = 1 / (2 * (n + lambda))

参数alpha控制Sigma点的散布范围，通常取1e-3到1之间；beta对高斯分布取2最优；kappa在状态维度高时取0或者3-n。电力系统状态维度往往不低，n大了之后(n+lambda)项会被拉大，Sigma点离中心过远，需要用更小的alpha来补偿。

3.3 什么时候UKF优势明显，什么时候并不值得

工程上的直觉判断是这样的：量测函数非线性强、状态维度在几十维以内时，UKF的精度和鲁棒性通常明显好过EKF，代价是计算量增大，因为每个Sigma点都要独立调用一次状态转移函数和量测函数。EKF只需要在单点算雅可比，计算量优势明显，但雅可比矩阵的推导、编程和验证是最费时间的环节。

如果系统接近线性、运行点变化不大，EKF的精度已经完全够用，没必要付出UKF几倍的计算开销。反过来，如果系统面临大扰动、故障或拓扑变化，EKF的一阶近似往往先撑不住，UKF因为做的是统计近似而不是局部线性化，适应性明显更强。在电力系统这种运行点经常大幅跳变的场景里，我倾向于优先把UKF作为基准滤波器，再用EKF作为计算量敏感场景下的替代方案。

4. Matlab实现：滤波器主循环、Sigma点传播与量测更新

4.1 初始化参数与数据准备

在写滤波主循环之前，先把状态维度、观测维度、噪声矩阵定义清楚。下面是一段我常用的初始化骨架，按二代模型写。

状态维度n等于发电机数量的两倍，观测维度m取决于量测配置。P0是初始协方差，代表对初值的不信任程度。如果初值来自上一个断面的静态状态估计结果，P0对角元素可以取小一些，比如1e-4；如果完全不知道初值，就取0.1以上，让滤波器前几步靠自己收敛。

% 基本参数 n = 2 * nG; % 状态维度：功角+转速 m = size(z, 1); % 量测维度 dt = 1 / 30; % 采样间隔，PMU典型30Hz % 初始状态与协方差 x0 = [delta0; omega0]; % 来自静态估计或潮流计算结果 P0 = 1e-4 * eye(n); % 噪声矩阵 Q = 1e-5 * eye(n); R = 1e-4 * eye(m); % 预分配状态轨迹 x_est = zeros(n, N); x_est(:, 1) = x0; P = P0;

Q和R的取值一定要结合物理量级。功角的方差和转速的方差不在一个量级，直接用相同数值会出问题。一般功角是弧度，偏差0.01弧度算很小；转速偏差2pi0.01弧度每秒也很小，但数值上这两者差了上百倍。Q的对角元素需要分别按功角、转速的实际波动范围来填。

4.2 EKF主循环代码骨架

以下是EKF核心循环的结构，实际使用时要替换状态转移函数f_sys、量测函数h_sys以及对应的雅可比矩阵函数。

for k = 2:N % ---------- 预测步骤 ---------- x_pred = f_sys(x_est(:, k-1), u(:, k-1), dt); A = F_jacobian(x_est(:, k-1), u(:, k-1), dt); P_pred = A * P * A.' + Q; % ---------- 更新步骤 ---------- H = H_jacobian(x_pred); S = H * P_pred * H.' + R; K = P_pred * H.' / S; innov = z(:, k) - h_sys(x_pred); x_est(:, k) = x_pred + K * innov; P = (eye(n) - K * H) * P_pred; % 保证协方差对称，避免数值误差累积 P = (P + P.') / 2; end

两个细节值得说明。雅可比矩阵H的计算点选在预测后的状态x_pred而不是上一时刻的估计值，这是工程惯例，因为量测更新是从预测点出发做校正的。矩阵求逆用P_pred * H.' / S而不是inv(S) * P_pred * H.'，因为除号在Matlab里调用的是数值稳定性更好的线性求解路径，显式求逆会引入额外的数值误差。

4.3 UKF主循环：Sigma点生成、传播与更新

UKF的主循环需要先生成Sigma点和权重，然后再进入递推过程。我直接给出一段可替换使用的Matlab代码骨架。

% Sigma点相关参数 alpha = 1e-3; beta = 2; kappa = 0; 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); Wm(2:end) = 1 / (2*(n + lambda)); Wc(2:end) = Wm(2:end); for k = 2:N % ---------- 预测步骤 ---------- sqrtP = chol((n + lambda) * P, 'lower'); X = zeros(n, 2*n+1); X(:, 1) = x_est(:, k-1); for i = 1:n X(:, i+1) = x_est(:, k-1) + sqrtP(:, i); X(:, i+n+1) = x_est(:, k-1) - sqrtP(:, i); end Y = zeros(n, 2*n+1); for i = 1:2*n+1 Y(:, i) = f_sys(X(:, i), u(:, k-1), dt); end x_pred = Y * Wm; P_pred = Q; for i = 1:2*n+1 d = Y(:, i) - x_pred; P_pred = P_pred + Wc(i) * (d * d.'); end % ---------- 量测更新步骤 ---------- Z = zeros(m, 2*n+1); for i = 1:2*n+1 Z(:, i) = h_sys(Y(:, i)); end z_pred = Z * Wm; Pzz = R; for i = 1:2*n+1 d = Z(:, i) - z_pred; Pzz = Pzz + Wc(i) * (d * d.'); end Pxz = zeros(n, m); for i = 1:2*n+1 dx = Y(:, i) - x_pred; dz = Z(:, i) - z_pred; Pxz = Pxz + Wc(i) * (dx * dz.'); end K = Pxz / Pzz; innov = z(:, k) - z_pred; x_est(:, k) = x_pred + K * innov; P = P_pred - K * Pzz * K.'; P = (P + P.') / 2; end

注意chol((n+lambda)*P, 'lower')要求输入矩阵严格正定。如果P因数值问题变成非对称或半正定，分解直接报错。前面EKF循环末尾的对称化处理，在UKF这里的作用就是为这一步兜底。如果仍遇到数值困难，可以在分解前做一个小扰动：P = P + 1e-12 * eye(n)，这在工程里很常见。

4.4 用数值雅可比校验手推公式

量测函数h_sys在EKF和UKF里必须完全一致，它是对比实验的公共基准。为了验证手推雅可比公式的正确性，推荐用中心差分方式计算数值雅可比，跟解析结果对比。

dx = 1e-6; H_num = zeros(m, n); for j = 1:n xp = x0; xm = x0; xp(j) = xp(j) + dx; xm(j) = xm(j) - dx; H_num(:, j) = (h_sys(xp) - h_sys(xm)) / (2*dx); end % 对比 H_num 和 H_analytic

如果某个元素的符号或者量级对不上，基本可以断定是手推公式有误。这个检查做一次，能省掉后面好几天的排查时间。

5. 算例对比：从追踪精度到鲁棒性，EKF和UKF的实际差异

5.1 算例设置与评价指标

以IEEE 39节点系统为例。这个系统俗称新英格兰系统，常用配置是10台发电机，对应状态维度20。如果想快速验证算法，可以先在3机9节点系统上做，状态维度只有6，趋势结论不变。

评价指标我一般用四组：状态估计的均方根误差、最大绝对误差、平均单步耗时、以及滤波器新息一致性。RMSE看整体水平，MAE看最坏情况，单步耗时看计算负担，新息一致性用NEES检验统计量来判断滤波器是否“过于自信”或“过于保守”。

仿真里有一点特别重要：所谓“真实状态”必须来自独立的仿真轨迹。先用一组真实的扰动模拟，比如负荷阶跃或短路故障，生成带噪声的量测数据，再把量测喂给滤波器。如果直接把滤波器自己的状态轨迹当真实值去算误差，那是在自己骗自己，结论没有意义。

5.2 典型精度与耗时对比

下面这张表给出的是这类对比中非常典型的结果模式，数值会随系统参数和量测噪声水平变化，但整体关系是稳定的。

UKF的精度优势在大扰动场景下更明显。功角RMSE从0.35度降到0.22度，看起来差距不大，但放到振荡阻尼分析和安全稳定校核里，0.1度的功角置信度差别会直接影响对系统稳定裕度的判断，这个是实打实的工程影响。

5.3 对粗差和突变场景的适应性差异

另一个值得关注的问题是粗差场景。PMU偶尔会给出远超噪声水平的错误量测。EKF的量测更新依赖雅可比矩阵在当前点的一阶等效，粗差出现时增益矩阵和残差的关系比较脆弱；UKF因为做的是统计拟合，对非高斯偏差有一定抗性。从实际仿真记录看，UKF在单点粗差后一般能在一到两个采样周期内把状态轨迹拉回来，EKF可能需要更长时间，严重时还会把后续状态预测带偏。

这里也要说句公道话：UKF并不能替代专门的坏数据检测与辨识模块。如果粗差持续存在，任何卡尔曼类滤波器都会被持续污染。工程实践里仍然建议在前面加一道数据合理性校验，比如量测突变幅值超过物理阈值就标记为可疑工况。

6. 调试经验：协方差整定、初值选择与代码效率优化

6.1 滤波器发散的第一排查顺序

滤波器发散是最常见的拦路虎。我的排查顺序很固定。

先画新息序列。如果新息均值明显不为零且持续有偏，多半是模型本身有错，不是参数问题。再看增益矩阵K的变化，如果K逐渐接近零并保持不变，说明滤波器已经完全信任模型、不再理会量测，这时任何量测错误都拉不回来。然后检查协方差矩阵是否正定，UKF里Cholesky分解失败是最直接的信号。最后看初值量级，功角初值差几十度也可能导致起始阶段发散，用上一个断面静态估计的结果做初值是最稳妥的策略。

6.2 Q和R的整定经验法则

Q和R的参数整定没有万能公式，但有一个实际可用的做法。先用一个比较宽松的R，也就是假设量测噪声偏大，跑一次仿真，记录残差或者预测误差的统计量，再反过来调整Q。也可以引入自适应思路：用滑动窗口内的新息协方差在线估计R。Matlab原型阶段我一般先固定Q和R，把算法验证通过之后，再考虑加入自适应模块。

还有一个容易被忽略的细节：把状态量和量测都转换到标幺值或者统一单位体系之后，Q和R的量纲一致，调参时直观很多。如果你用弧度表示功角、用标幺值表示功率、用有名值表示电压，这三个量级的数值差异会让协方差矩阵变成病态，实际效果就是参数怎么调都不对。

6.3 从Matlab原型到实际部署的注意点

Matlab里跑通只是第一步。真实工程中的动态状态估计模块往往有实时性要求，下面几个点值得注意。

数组预分配。主循环涉及的历史轨迹存储，一定要提前分配好空间，不要靠在循环里动态扩展数组，那在长时间运行时会非常慢。Sigma点传播很难完全向量化，因为每个Sigma点都要独立调用一次非线性函数，但可以把函数入口的参数结构提前算好，减少重复计算。如果机器是多核，可以对Sigma点传播使用并行循环，但要注意并行计算对函数作用域和随机依赖的限制。最后，MATLAB Coder可以把滤波主循环转成C代码，但对动态内存分配和匿名函数限制较多，所以早期写代码时就要有意识地避免在循环内部使用匿名函数和大小可变的数组。

7. 个人实验心得：EKF和UKF怎么选

我在实际做这类对比实验时的一个核心体会是：EKF和UKF在电力系统动态状态估计里没有绝对的谁替代谁，更多的是看具体场景。

如果系统只是常规负荷波动，运行点变化不大，EKF的计算效率优势值得保留，而且实现代码更短，调试成本低。如果系统经常面临大扰动、量测配置复杂、可能需要在线处理粗差，UKF的鲁棒性和精度优势更明显，多出来的计算量在现在的主流硬件上通常可以接受。

有一个小技巧值得分享：不论用哪种滤波器，把状态转移函数和量测函数的输入输出维度做成可配置的，在同一个算例上先跑EKF再跑UKF，对比两条新息序列曲线。新息均值、波动范围的一致性，能帮你快速定位模型实现里的隐性错误，这比单纯看RMSE数字直观得多。我自己的习惯是把这条新息曲线作为每个新算例的“体检项”，先看它再谈其他指标。