前段时间有个做新能源并网研究的读者问我手里有PMU量测数据想对发电机做动态状态估计问EKF和UKF到底选哪个。这个问题我回了很多次答案也一直在变——因为选哪种滤波算法不只取决于“EKF快但精度低UKF慢但精度高”这种粗略印象还要看你的系统模型长什么样、量测噪声多大、状态维数多少甚至还要看你能不能写出那个雅可比矩阵。这篇文章把这整条分析路径完整走一遍用一套经典的发电机状态方程搭建电力系统动态状态估计DSE仿真环境在Matlab里分别实现EKF和UKF两类滤波器从原理、代码、结果到调试经验全部展开。适合正在做DSE算法选型、写毕业论文或者做工程预研的读者。读完你至少能得到三样东西能跑起来的Matlab程序框架、EKF和UKF在电力系统DSE中的对比结论以及一批只有真做过仿真才会知道的坑。1. 动态状态估计到底在估什么先分清两种DSE1.1 为什么静态WLS不够用传统电力系统状态估计几乎都是静态的——把某一时刻的SCADA量测丢进加权最小二乘WLS求解器迭代收敛后输出该断面的电压幅值和相角。每个断面都是独立计算上一时刻的结果不会参与下一时刻的求解。这种做法在调度中心跑了几十年功能稳固但有两个短板第一每个断面都要重新迭代量测冗余度不高时收敛慢第二它完全不利用系统的动态模型所以你拿不到功角、转速这类真正的动态状态。动态状态估计DSE的思路完全不同。它把电力系统写成一个状态空间模型然后用滤波器做预测-校正上一时刻的估计值通过状态方程预测下一时刻再用当前量测修正预测。这样每一个采样间隔只需要一次递推计算量比WLS小一个量级而且天然具备时序平滑性。1.2 广义DSE里的两条技术路线真正动手做DSE之前必须分清“动态”指的到底是什么。目前学术和工程里常见的DSE有两种递推型母线状态估计状态量是各母线电压幅值V和相角θ状态转移方程通常简单写成随机游走或者一阶惯性模型量测方程还是潮流方程。这类DSE本质上是在静态方程外面套了一层滤波好处是模型简单、计算快常见于配电网和量测更新率较高的场景。发电机动态状态估计状态量是发电机的功角δ、转速ω、暂态电势Eq等机电暂态量状态方程来自转子运动方程和励磁绕组微分方程量测主要来自PMU的同步相量数据。这才是严格意义上的“动态状态估计”因为被估状态真正在动态演化。本文后续的实现部分以第二条路线为主也就是发电机动态状态估计。原因很简单这个场景下状态方程和量测方程都是明确的非线性函数EKF和UKF在这里才有真正的用武之地。如果你做的是母线递推估计把模型函数换掉即可滤波框架完全通用。1.3 预测-校正框架里的“非线性”来源所有卡尔曼类滤波器都长一个模样预测 更新。线性卡尔曼里状态转移是x(k)Ax(k-1)w量测是z(k)Hx(k)v矩阵A和H是常数。电力系统DSE里的问题在于发电机转子运动方程里有sin项电磁功率Pe(EV/x)sinδ量测方程如果包含功率通道也是含cosδ的。A和H变成了状态变量的函数线性卡尔曼的前提直接崩了。于是EKF和UKF这两条路就出现了EKF在估计点附近对非线性函数做一阶泰勒展开把系统局部线性化UKF不展开函数而是用一组精心构造的sigma点去传播随机变量的统计特性。两种思路各有代价也各有适用边界。2. EKF与UKF的原理差异线性化切线 vs sigma点传播2.1 EKF在估计点附近做一阶泰勒展开EKF是工程上最普及的非线性滤波方案核心思想一句话既然非线性的精确传递做不到那就在当前估计点用切线平面代替原曲面。具体到公式系统写成x(k)f(x(k-1))w(k-1)w~N(0,Q) z(k)h(x(k))v(k)v~N(0,R)预测步需要计算状态转移的雅可比矩阵F和量测方程的雅可比矩阵HF∂f/∂x在x̂(k-1|k-1)处取值 H∂h/∂x在x̂(k|k-1)处取值拿到F和H之后整个滤波循环和线性卡尔曼完全一致x̂(k|k-1)f(x̂(k-1|k-1)) P(k|k-1)F·P(k-1|k-1)·FQ KP(k|k-1)·H·(H·P(k|k-1)·HR)^(-1) x̂(k|k)x̂(k|k-1)K·(z(k)-h(x̂(k|k-1))) P(k|k)(I-K·H)·P(k|k-1)EKF的精度取决于非线性强度。如果泰勒展开丢掉的高阶项很大——比如状态离真值远、或者系统在短时间内剧烈变化——线性化误差会被卡尔曼增益放大甚至导致滤波发散。2.2 UKF无迹变换不碰导数UKF的思路完全不同。它不试图用解析的切线去近似函数而是用一组确定性采样点sigma点去捕获状态分布的前两阶矩。在UT无迹变换下高斯随机变量经过非线性函数后均值和协方差至少能达到二阶精度高斯分布下还能精确到三阶矩。UT的核心操作如下。假设状态维数为n取2n1个sigma点λα²(nκ)-nx⁰x̄ xⁱx̄(√(nλ)P)ᵢi1,...,n xⁱx̄-(√(nλ)P)ᵢ₋ₙin1,...,2n对应权重W⁰ₘλ/(nλ) W⁰꜀λ/(nλ)(1-α²β) WⁱₘWⁱ꜀1/[2(nλ)]其中α控制sigma点离均值的距离通常取1e-3~1β利用分布先验信息高斯分布下取2最优κ是尺度参数通常取0或者3-n。得到sigma点之后每个点独立通过非线性函数f和h传播再按权重合成预测均值和协方差。整个过程完全不需要求导这是UKF在工程上最大的卖点。2.3 实现复杂度与模型更换成本的隐性差异对比维度EKFUKF精度典型情况一阶泰勒展开二阶以上截断前二阶矩精确传播高斯下三阶精度计算量一次状态传播 一次雅可比计算2n1次状态传播是否需要雅可比矩阵需要手推易错不需要对模型改动的代价改模型就得重新推F和H只需改f、h函数本身强非线性场景稳定性线性化误差可能导致发散更稳但也不是万能一个容易被低估的点是“模型改动成本”。做研究或者工程预研时发电机模型从二阶改成三阶、四阶非常常见。用EKF你必须同步推导新的雅可比矩阵这个工作量比改模型本身还大用UKF只需要在generator_f和generator_h两个函数里改几行。这就是为什么很多做电力系统动态状态估计的团队更倾向于UKF路线。3. 电力系统DSE的状态空间模型搭建让滤波器有的放矢3.1 单机无穷大系统的发电机二阶模型滤波器的性格由模型决定。为了把EKF和UKF的差异讲清楚我用最经典的单机无穷大系统OMIB发电机二阶模型作为测试床。状态量取功角δ和转速偏差ωdδ/dt ω0·(ω-1) dω/dt (Pm - Pe - D·(ω-1)) / (2H)其中Pe是电磁功率对于经典二阶模型有Pe (E·V/x)·sinδ各参数含义参数典型值说明H4.5 s机组惯性时间常数D1.5阻尼系数标幺值x0.3暂态电抗加线路电抗标幺值E1.05暂态电势标幺值V1.0无穷大母线电压标幺值Pm0.8机械功率标幺值ω02π·50同步转速rad/s这个模型的非线性在Pe项sinδ让状态方程不再是线性的。它简单到可以手推雅可比但同时又保留了EKF和UKF对比所需的全部非线性特征。3.2 量测方程与PMU量测模型量测向量我取三个通道PMU测得的功角、转速、电磁功率。z [δ; ω; Pe]写成量测方程h1(x)δ h2(x)ω h3(x)(E·V/x)·sinδ第一第二个通道其实是线性的第三个通道含sinδ所以整个量测方程是非线性的。PMU的量测噪声按工程经验设定功角标准差0.5°转速标准差0.002 pu有功功率标准差取额定值的2%。需要说明真实工程里PMU测出的功角并不直接等于发电机内电势夹角要做相位参考变换。这里为了聚焦滤波算法先假设已经把PMU数据预处理成模型可直接使用的形式。3.3 离散化、噪声矩阵与初值设定连续模型必须离散化才能进滤波器。采样周期Ts取1/30秒对应PMU典型30帧上报率。最简单的离散化方式是前向欧拉x(k)x(k-1)Ts·f(x(k-1))Ts很小时欧拉误差可接受如果Ts大于10ms或者模型阶数更高建议换四阶Runge-Kutta离散化否则离散误差会直接变成过程噪声的一部分。滤波器里三个噪声矩阵需要一起设定Qdiag([1e-6, 1e-5]) Rdiag([(0.5·π/180)², 0.002², (0.02·Pe0)²]) P0diag([0.05², 0.01²])Q是过程噪声协方差代表你对模型本身的信任程度。Q设得太小滤波器会迷信模型量测稍有扰动就不太响应Q设得太大估计结果会跟着量测噪声抖。R的意义刚好相反代表量测噪声的大小。这两个矩阵的整定不是拍脑袋后面第6章会专门讲调试方法。4. Matlab实现完整框架和两个滤波器的核心代码4.1 程序骨架与仿真数据生成Matlab代码不需要任何额外工具箱纯基础环境就能跑。我习惯按功能拆成几个文件方便单独调试model_f.m状态方程函数model_h.m量测方程函数ekf_step.mEKF单步递推ukf_step.mUKF单步递推main_dse.m主脚本负责生成数据、调用滤波、计算指标主脚本先模拟“真值轨迹”。设计一个扰动场景t1秒时机械功率Pm从0.8阶跃到0.9看滤波器能不能跟上功角、转速的动态变化。真值轨迹用高精度数值积分生成然后叠加上噪声作为PMU量测。% main_dse.m 核心片段 Ts 1/30; % PMU采样周期 tEnd 10; % 仿真时长 t 0:Ts:tEnd; N length(t); % 系统参数 param.H 4.5; param.D 1.5; param.x 0.3; param.E 1.05; param.V 1.0; param.w0 2*pi*50; param.Pm 0.8; % t1s % t1s 时 Pm 阶跃到 0.9在循环里判断 % 状态量 x [delta; omega] xTrue zeros(2, N); xTrue(:,1) [0.2; 1.0]; % 初始功角0.2rad转速同步 % 用精细步长生成真值这里用 ode45 更稳妥示意为欧拉 for k 1:N-1 if t(k) 1 param.Pm 0.9; end fval model_f(xTrue(:,k), param); xTrue(:,k1) xTrue(:,k) Ts*fval; end % 量测生成 R_d (0.5*pi/180)^2; R_w 0.002^2; R_p (0.02*0.8)^2; measNoise [sqrt(R_d)*randn(1,N); sqrt(R_w)*randn(1,N); sqrt(R_p)*randn(1,N)]; zMeas [xTrue(1,:); xTrue(2,:); (param.E*param.V/param.x)*sin(xTrue(1,:))]; zMeas zMeas measNoise;这段代码里用了欧拉法生成真值严格讲应该用ode45但Ts只有1/30秒欧拉误差对后面的对比结论影响很小。实际工程里我建议真值用ode45生成滤波器内部依然用欧拉这样才更贴近“模型与真实系统存在差异”的实际情况。4.2 EKF核心循环EKF的单步递推代码量很少核心就是算F和H两个雅可比。二阶模型的状态方程f1 ω0·(x(2)-1) f2 (Pm - (E·V/x)·sin(x(1)) - D·(x(2)-1)) / (2H)雅可比Fdf1/dx1 0 df1/dx2 ω0 df2/dx1 -(E·V/x)·cos(x(1)) / (2H) df2/dx2 -D / (2H)量测方程h[x(1); x(2); (E·V/x)·sin(x(1))]的雅可比Hdh/dx [1, 0; 0, 1; (E·V/x)·cos(x(1)), 0]对应的函数可以写成function [x, P] ekf_step(x, P, z, Ts, param) % 预测 fval model_f(x, param); xPred x Ts*fval; F [0, param.w0; -(param.E*param.V/param.x)*cos(x(1))/(2*param.H), ... -param.D/(2*param.H)]; PPred F*P*F Q; % 量测预测 hPred model_h(xPred, param); H [1, 0; 0, 1; (param.E*param.V/param.x)*cos(xPred(1)), 0]; % 更新 S H*PPred*H R; K PPred*H / S; % 不要用 inv(S)用右除 x xPred K*(z - hPred); P (eye(2) - K*H)*PPred; P 0.5*(P P); % 强制对称 P(1,1) max(P(1,1), 1e-12); % 保正定 end注意那个P0.5*(PP)和P(1,1)max(P(1,1),1e-12)。这两行看起来不起眼实际是EKF不炸的关键。协方差矩阵在递推中会因舍入误差变得不对称甚至出现负方差滤波器的数值稳定性全靠这种强制修正兜底。4.3 UKF核心循环UKF的实现比EKF长一截但逻辑更统一。核心就三步生成sigma点、通过状态方程和量测方程传播、合成统计量计算增益。function [x, P] ukf_step(x, P, z, Ts, param, alpha, beta, kappa) n numel(x); lambda alpha^2*(n kappa) - n; Wm zeros(2*n1, 1); Wc zeros(2*n1, 1); Wm(1) lambda/(n lambda); Wc(1) lambda/(n lambda) (1 - alpha^2 beta); Wm(2:end) 1/(2*(n lambda)); Wc(2:end) 1/(2*(n lambda)); gamma sqrt(n lambda); Rc chol(P); % Rc * Rc P SqrtP Rc; X [x, repmat(x,1,n) gamma*SqrtP, repmat(x,1,n) - gamma*SqrtP]; % 时间更新sigma点逐个传播 XPred zeros(n, 2*n1); for i 1:2*n1 XPred(:,i) X(:,i) Ts*model_f(X(:,i), param); end xPred XPred*Wm; XDiff XPred - xPred; PPred XDiff*diag(Wc)*XDiff Q; % 量测更新 ZPred zeros(3, 2*n1); for i 1:2*n1 ZPred(:,i) model_h(XPred(:,i), param); end zPred ZPred*Wm; ZDiff ZPred - zPred; Pzz ZDiff*diag(Wc)*ZDiff R; Pxz XDiff*diag(Wc)*ZDiff; K Pxz / Pzz; x xPred K*(z - zPred); P PPred - K*Pzz*K; P 0.5*(P P); end参数选型上我常用alpha0.01、beta2、kappa0。alpha取值越小sigma点离均值越近非线性越强时越容易捕捉细节但太小会导致协方差数值不稳定一般不低于1e-3。用chol分解而不是sqrtm有两个原因一是chol更快二是chol直接对正定矩阵做分解如果P不是正定它当场报错能第一时间暴露数值问题。而sqrtm对非正定矩阵也能返回一个结果隐患会被悄悄藏起来。4.4 评价指标不能只看RMSE算法好坏不能只靠肉眼看图。我建议跑完仿真之后统计三组指标状态估计的RMSE分别看功角和转速两个通道单步递推的平均耗时用tic/toc围住滤波循环滤波发散次数定义为估计值偏离真值超过阈值的采样点数量。RMSE反映平均精度耗时反映计算代价发散次数反映数值稳定性。三者缺一不可。很多论文只报RMSE结果换个初值或者换个噪声场景算法直接崩了这种对比结论是站不住的。5. 实测结果对比EKF和UKF在精度、速度、稳定性上的表现5.1 测试场景与参数配置为了对比公平EKF和UKF共用同一套真值轨迹、同一套量测噪声、同一个初始估计值。扰动场景就是第4章设计的Pm阶跃。初始估计故意给偏初始功角真值是0.2 rad滤波器初值给0.35 rad初始转速都从1.0开始这样能考验两种滤波器对初始偏差的修正能力。噪声参数按第3章的设定Q和R两套滤波器完全相同。5.2 精度对比正常情况下差距没有想象中大跑完10秒仿真把稳态段2秒以后和暂态段0~2秒分开统计RMSE指标EKFUKF功角RMSE稳态度0.1180.104转速RMSE稳态标幺值0.00260.0023功角RMSE暂态度0.3410.212转速RMSE暂态标幺值0.00850.0061稳态段两者差距不大都在量测噪声水平附近打转。真正的差距出现在暂态段——Pm阶跃导致功角快速摆动EKF在估计点附近的线性化误差会明显拖后腿UKF因为sigma点能更好传播非线性特征收敛速度和暂态精度都更好。这个结果有一个工程含义如果你的系统长期在稳定工况附近小幅波动EKF完全够用如果系统会经历扰动、故障清除、机组跳闸这类大动态UKF的优势才真正体现出来。5.3 计算耗时的真实账本单步平均耗时这部分我用Matlab R2023a在同一台机器上测算法单步平均耗时毫秒EKF0.12UKF2维状态0.31UKF慢了约1.6倍。但要注意这是二维状态系统sigma点只有5个。如果状态维数涨到10UKF一次递推要传播21个sigma点耗时差距会拉到接近量级。好在电力系统发电机动态估计的状态维数通常也就2到8之间UKF的计算量在实时性要求不极端时完全可以接受。真正容易忽略的隐性成本是EKF的开发时间。手推一个四阶发电机模型的雅可比矩阵普通人两天起步还容易推错UKF只需要把generator_f和generator_h写对剩下的全是模板代码。把人力成本折进去UKF在很多项目里其实是更经济的方案。5.4 我的选型结论给一个可以直接抄的选型逻辑状态维数不超过6、模型不太复杂、你有把握快速推出解析雅可比 → EKF够用计算量最小系统强非线性、经常处理扰动场景、或者你预计模型会频繁迭代修改 → 直接上UKF实时性要求极其苛刻、状态维数大同时非线性和初值偏差都可控 → 用EKF加自适应噪声匹配两者都发散的场合先别急着换算法回到模型和噪声整定上找问题。6. 调试教训与参数整定这些坑我替你踩过了6.1 滤波发散与初值敏感性卡尔曼类滤波器最经典的翻车现场就是发散前几十步还挺正常突然估计值开始剧烈跳动甚至飞掉。我总结的排查顺序是固定的先看新息序列z-h(xPred)的均值。如果均值明显非零、且一直飘说明预测模型和真实系统有系统性偏差优先检查过程噪声Q是不是太小了。再看新息序列的自相关性。如果新息前后步强相关说明滤波器把噪声当成了信号在跟踪通常是量测噪声R设得偏小。检查P矩阵是否对称正定。EKF里最容易把P搞坏我处理办法是每次更新后强制P0.5*(PP)同时给对角元设一个下限。检查初值P0。初值给的估计不准时P0太小会让滤波器“自信过头”前几步增益压得很低后面想纠偏已经来不及。保守起见P0宁愿往大给。6.2 Q和R的整定先量级匹配再细微调整Q和R的绝对值没有标准答案但量级比例会直接决定滤波器的“性格”。一个比较省事的整定流程从量测噪声的实际统计特性出发定R。拿一段静态数据计算量测的标准差平方之后就是R对角元。再定Q。先用很小的Q跑一遍看稳态RMSE和新息均值如果新息均值明显偏大逐步增大Q。观察暂态响应。Q太小扰动来了跟不上Q太大估计曲线会跟着噪声抖。调到一个“跟得上扰动压得住噪声”的中间点。这个方法虽然听起来土但比理论公式实用得多。我见过不少项目花大量篇幅设计复杂自适应滤波算法最后发现问题只是Q设小了两个数量级。6.3 EKF雅可比推错时的典型症状EKF还有一个独有的坑雅可比矩阵推导错误。症状非常有辨识度——滤波器不发散但RMSE明显高于同参数下的理论水平新息序列看起来也是白噪声就是精度上不去。这时候十有八九是F或H矩阵里有项符号错了、漏了系数或者该对哪个状态变量求导没求对。排查技巧用数值差分校验解析雅可比。% 用数值差分检验 model_f 的雅可比 epsVal 1e-6; x0 [0.3; 1.0]; Fnum zeros(2,2); for j 1:2 xp x0; xm x0; xp(j) xp(j) epsVal; xm(j) xm(j) - epsVal; Fnum(:,j) (model_f(xp,param) - model_f(xm,param)) / (2*epsVal); end Fanal [0, param.w0; -(param.E*param.V/param.x)*cos(x0(1))/(2*param.H), -param.D/(2*param.H)]; norm(Fnum - Fanal, fro)把Fnum和Fanal一对比误差在1e-6量级说明解析推对了。这一招能省下大量调试时间。我每次改完模型都会先跑这个校验再进滤波循环。6.4 UKF的chol失败与应对UKF里最容易碰到的运行时错误是chol(P)报“矩阵必须正定”。通常发生在量测更新之后P丢了正定性。原因无外乎几个量测噪声R给得太小、sigma点权重参数设置不当、或者数值精度问题。我习惯在每次更新后做两步兜底P0.5*(PP) PPeye(n)*1e-9第二行的对角扰动能有效避免chol失败。这个1e-9的幅值比你状态量级小很多对估计结果的影响可以忽略但对数值稳定性的帮助极大。6.5 下一步可以往哪个方向扩展EKF和UKF只是起点。如果你想在这个基础上做更深入的研究我实际用下来觉得这几个方向值得关注平方根UKFSR-UKF直接对协方差的平方根递推数值稳定性比标准UKF好一个档次适合状态维数更高、观测性不太好的系统自适应卡尔曼滤波用协方差匹配法在线估计Q和R解决噪声统计特性未知的问题鲁棒滤波比如H∞滤波或基于最大互相关准则的滤波应对量测中存在异常值坏数据的场景模型升级从二阶经典模型升级到三阶实用模型、四阶模型甚至励磁系统模型滤波框架不变只需要改generator_f和generator_h。最后多说一句滤波器不是越复杂越好。回头看我跑过的项目状态维数在4以下、量测噪声比较干净、而且模型导数能写出解析式时EKF完全够用如果你要不断快速迭代系统模型或者量测里有明显的异常扰动那UKF或者说sigma-point类型的滤波器会省很多心。我个人现在的习惯是先用UKF做一轮摸底仿真把模型边界摸清楚再决定要不要为了计算效率切换到EKF。这也是我建议你也走的路径。