企业官网建设流程全解析

卡尔曼滤波家族在电力系统动态状态估计里，一直是讨论热度比较高的方向。尤其是扩展卡尔曼滤波（EKF）和无迹卡尔曼滤波（UKF）这两兄弟，一个靠线性化硬啃非线性，一个靠Sigma点采样“曲线救国”，很多做电力系统仿真和实时监控的朋友都在这两者之间反复权衡过。我最早接触这块是在做单机无穷大系统的功角估计时，手里有一堆Matlab脚本，但真正把EKF和UKF从公式变成能跑、能对比、能出图的代码，还是踩了不少坑。这篇文章就把我完整跑通的实现过程、原理拆解和调参经验整理出来，希望能帮到正在做电力系统动态状态估计课题或者毕业设计的同学。

1. 为什么电力系统需要“动态”状态估计：从静态断面到实时跟踪

很多人刚接触状态估计时，脑子里默认是电力系统里那个经典问题：给定SCADA量测，用加权最小二乘（WLS）去算节点电压幅值和相角。这就是所谓的静态状态估计。静态估计的隐含假设是系统运行在某个稳态断面，量测一帧一帧地独立处理，每一帧之间没有时间上的耦合关系。在早期电网结构相对简单、负荷平稳的场景下，这个假设是够用的。

但现在的系统早就不一样了。新能源大规模接入，风机和光伏出力波动剧烈；直流输电投运后，功率传输走廊的动态行为更复杂；系统的等效惯量在下降，功角振荡、电压失稳这类动态过程出现的频率和烈度都在上升。在这种背景下，“下一个时刻的状态”不再是独立随机变化的，它严格由当前状态和系统动力学方程决定。这时候如果还用静态估计去逐帧处理，本质上是在扔掉一个极其重要的信息——系统动态模型本身。动态状态估计的核心价值，就是把这个“状态随时间演化”的规律纳入滤波框架，让估计器不仅能看懂当前断面，还能利用上一时刻的信息对下一时刻做预测，再用量测修正预测，从而实时跟踪系统的动态轨迹。

具体到物理对象，动态状态估计主要关注发电机的动态过程：功角δ、角频率偏差Δω、q轴暂态电动势E'q，这些量共同刻画了发电机在扰动后的摇摆行为。对这些状态做实时估计，意义不只是“看着好看”，而是为了给后续的紧急控制、切机切负荷、阻尼控制提供可靠的实时反馈信号。很多控制策略在实施时，如果状态反馈信号本身是滞后或者含噪的，控制效果会大打折扣。

从信息融合的角度看，动态状态估计把两类信息拧在了一起：一是系统模型给出的“预测”，二是PMU/SCADA量测给出的“修正”。预测精不精、修正强不强，直接决定了估计器的性能。EKF和UKF就是两种实现这种预测-修正循环的滤波策略，区别只在于“如何把非线性模型塞进卡尔曼框架”的手法不同。

2. 动态状态估计的建模基础：状态方程与量测方程到底怎么列

在进入EKF和UKF算法细节之前，先把模型立起来。我用的是最经典的单机无穷大系统（SMIB），这几乎是动态状态估计研究的“Hello World”。虽然简单，但它的非线性特征、状态耦合、量测映射关系和多机系统在数学结构上完全一致，跑通它再扩到多机就是体力活。

2.1 状态向量的选取：发电机经典三阶模型

发电机动态模型可以取到不同阶数。工程上常用的是三阶模型，状态向量取：

x = [δ, Δω, E'q]

三个状态的物理含义很直观：

• δ是发电机的转子功角，单位是rad；
• Δω是转子角速度相对于同步转速的偏差，单位是rad/s；
• E'q是q轴暂态电动势，单位是p.u.。

这三者构成了一个完整的“机电暂态+励磁动态”的最小集合。如果你在文献里看到把汽轮机调速器状态也塞进来的，那就是更高阶模型，原理完全一样，只是状态维数多了，滤波器的计算量会涨。

2.2 连续时间状态方程的离散化处理

发电机三阶模型的连续时间微分方程如下：

dδ/dt = ωb * Δω dΔω/dt = (Pm - Pe - D*Δω) / (2H) dE'q/dt = (Efd - E'q - (xd - x'd) * Id) / T'd0

这里的参数里，ωb是同步角速度基准值（一般是2π50或者2π60），Pm是机械功率，Pe是电磁功率，D是阻尼系数，H是惯性时间常数，Efd是励磁电压（恒定励磁时为一个常量），xd和x'd分别是同步电抗和暂态电抗，T'd0是励磁绕组开路时间常数，Id是d轴电流。

量测方程建立在电气量上。常见的量测包括机端电压幅值Vt、有功功率Pe、无功功率Qe，它们和状态的关系是：

Pe = (E'q * V∞ / x'dΣ) * sin(δ) Qe = (E'q^2 / x'dΣ) - (E'q * V∞ / x'dΣ) * cos(δ) Vt的表达式稍微复杂一点，和E'q、δ、系统等值电抗都有关。

可以看到，量测方程相对于状态是高度非线性的，这正是EKF和UKF发挥价值的地方。如果模型是线性的，直接用标准卡尔曼滤波就行，也就不需要在这篇文章里讨论这三兄弟的区别了。

在Matlab里搭建模型时，我的做法是先把连续方程写成函数句柄，再用欧拉法离散化。欧拉法在采样周期足够小（比如0.01s以下）的时候精度足够，而且实现最简单。如果你的采样时间较大（比如0.02s以上），建议改成四阶龙格库塔（RK4）来离散状态转移，否则预测误差会偏大，滤波精度会受拖累。用RK4会多几次函数求值，但动态估计的采样频率通常就是几十到上百赫兹，这个计算开销完全可接受。

3. EKF在电力系统动态估计里的实现：线性化这一步做对了，后面才有意义

EKF的基本思想一句话就能说清楚：既然卡尔曼滤波要求线性高斯系统，那我就在每个时刻把非线性函数在估计点附近做一阶泰勒展开，得到一个“瞬时线性化”的模型，然后套标准卡尔曼流程。

3.1 雅可比矩阵的两种求法：解析推导和数值差分

EKF在实现时最麻烦的就是算两个雅可比矩阵：状态转移矩阵F和量测矩阵H。很多教材上只写“求雅可比”，真上手时会发现这一行字能让人折腾一天。

F矩阵是状态函数f对状态x的偏导，H矩阵是量测函数h对状态x的偏导。对三阶发电机模型来说，解析推导虽然繁琐但可行。比如对Pe求δ的偏导，就是dPe/dδ = (E'q*V∞/x'dΣ)*cos(δ)，这个写起来不算太难。但一旦状态维数涨到10维以上，解析推导雅可比矩阵就变成了非常痛苦的差事，而且容易出错。

我个人的建议是：在验证阶段用解析雅可比，确认算法跑通之后，再换成数值差分，代码更简洁，扩展性更好。数值差分在Matlab里实现很直接，可以用中心差分：

for i = 1:n px = x0; px(i) = x0(i) + h; mx = x0; mx(i) = x0(i) - h; F(:,i) = (f(px) - f(mx)) / (2*h); end

这里的步长h需要谨慎选择，过大会导致截断误差，过小会引入舍入误差。我的经验值是取 sqrt(eps) * max(1, abs(x0(i)))，效果比较稳定。

3.2 EKF预测步：状态与协方差的时间更新

预测步的流程在数学上很干净：用当前时刻的后验估计值x_hat(k)代入非线性状态转移函数，得到下一时刻的先验估计；误差协方差则通过雅可比矩阵线性传播：

x_pri(k+1) = f(x_est(k)) P_pri(k+1) = F * P_est(k) * F' + Q

这里Q是过程噪声协方差矩阵。有人会问，为什么协方差的预测要经过F矩阵左乘右乘？打个不严谨但好理解的比方：如果你对一个向量做线性变换，那这个向量的“不确定范围”（近似看作一个椭球）也会被同一个变换拉伸、旋转。F矩阵就是那个变换的局部代言人，所以协方差要两边同时乘F。

3.3 EKF更新步：卡尔曼增益与后验修正

更新步骤完全是标准卡尔曼滤波的形式：

K = P_pri * H' * (H * P_pri * H' + R)^(-1) x_est(k+1) = x_pri + K * (z - h(x_pri)) P_est(k+1) = (I - K * H) * P_pri

在实际代码里，计算增益时千万不要显式求逆。用Matlab的右除运算符会更稳定：K = (P_pri * H') / (H * P_pri * H' + R)。这个细节在P矩阵病态时能救你一命。

更新步的直觉也值得说清楚：Kalman增益K本质上是一个加权系数，它权衡“模型的预测”和“量测的修正”谁更可信。如果R特别小（量测很准），K会偏向量测；如果Q特别大（模型不可信），同样是K偏向量测。这个平衡关系在调参阶段的理解至关重要。

3.4 EKF的实际局限：强非线性场景下的“硬伤”

用EKF在电力系统里做动态估计，最容易翻车的场景是大扰动后的暂态过程。功角在故障清除后会发生大范围摆动，此时线性化点的邻域很小，一阶泰勒展开很快就撑不住了，误差可能被急剧放大。我做过一个实验：在三相短路故障场景下，如果故障持续时间超过0.1s，EKF的功角估计误差有时会冲到1度以上，这在严苛的动态安全分析里是不可接受的。EKF计算量小、实现灵活，但它把宝押在“局部线性化足够准”上，这个假设在强非线性时非常脆弱。

4. UKF实现细节：用Sigma点“无迹变换”绕开雅可比矩阵

UKF的出发点完全换了个思路：我不在估计点附近去做线性化，而是直接找一组精心挑选的采样点（Sigma点），把这组点分别通过非线性函数传播，再根据传播后的点集重新统计均值与协方差。这个操作被称为无迹变换（Unscented Transform, UT）。

4.1 为什么Sigma点比随机采样更优雅

你可能第一时间想到蒙特卡洛：扔几万个随机点传播过去，均值协方差也能统计出来吧？理论上可以，但粒子滤波的教训告诉我们，随机采样的计算量和估计精度之间的权衡很尴尬。UT的精致之处在于，Sigma点是确定性选取的，2n+1个点就能抓住高斯分布的一阶矩和二阶矩信息，传播后得到的均值和协方差可以精确到非线性函数的二阶泰勒展开项，这一精度已经超过EKF的一阶线性化，而计算量只多了几次函数求值。

4.2 Sigma点生成与权重计算

对n维状态（我这里n=3），选取2n+1 = 7个Sigma点。比例对称采样法的经典做法如下：

计算矩阵平方根：S = sqrt((n+λ) * P)，其中λ = α^2*(n+κ) - n。然后用chol()函数分解协方差矩阵。

Sigma点这样生成：

• 第0个点就是当前状态均值；
• 第1到第n个点等于均值加上S的第i列；
• 第n+1到第2n个点等于均值减去S的第i列。

权重计算要分两套：一套用于求均值，一套用于求协方差。均值权重里第0个点是 λ/(n+λ)，协方差权重第0个点是 λ/(n+λ) + (1-α^2+β)，其余点的两组权重都是 1/(2(n+λ))。这里的经验参数取值：α控制Sigma点离均值的距离，建议取1e-3到1之间，系统越非线性取越小；β对高斯分布取2最优；κ一般取0，保证半正定。

很多文献会说这个步骤“简单”，但真实写代码时有一个很隐形的坑：如果P矩阵在迭代过程中因为数值误差变得不再正定，chol()会直接报错。防这个问题的办法有两个，一是给P加一个微小对角阵（比如1e-12*I），二是直接用sqrtm()代替chol()，虽然慢一点，但能兼容半正定矩阵。

4.3 UKF的状态预测与量测更新

UKF的预测步很直白：把7个Sigma点逐一丢进状态转移函数，得到一组传播后的点，然后按照权重加权平均，得到先验状态估计；对每个点减去先验均值的加权外积求和，再加上过程噪声Q，就得到先验协方差。

量测更新的关键选择是：用“重组Sigma点”还是“复用预测Sigma点”？严格的做法是拿先验状态重新生成Sigma点，再丢进量测函数；但更省计算量的做法是直接把预测步算出来的传播后状态点丢进量测函数。两种做法在理论上都有文献支持，我实际测试的结果是差异极小，为了代码简洁通常选择复用。

量测均值、量测协方差和状态-量测互协方差算出来后，计算UKF增益：

K = P_xz / P_zz x_est = x_pri + K * (z - z_mean) P_est = P_pri - K * P_zz * K'

这里P_xz是状态与量测的互协方差，P_zz是量测协方差。整个流程里完全没有求导，这就是UKF比EKF最大的工程优势——你把状态方程和量测方程变成复杂的查表函数、甚至Simulink仿真模型，UKF也能照样滤波。

4.4 UKF在不同场景下的稳定性表现

我跑过两个典型场景对比：平缓负荷波动和三相短路故障。在平缓场景里，EKF和UKF的估计误差都在0.01度量级，差异小到几乎看不出区别；但在短路故障后的功角大幅振荡场景，UKF的跟踪能力明显更稳，功角最大估计误差比EKF小一个数量级。代价就是计算量：每个时刻要执行至少7次状态函数求值和7次量测函数求值，比EKF的几次雅可比计算+两次函数求值要贵。但电力系统动态估计的采样频率通常是工频的1到2个周波一次，也就是10ms到20ms一次，这点计算量在现在的工控机和工作站上完全不是瓶颈。

5. Matlab仿真实验：从真值生成到RMSE对比的完整流程

这一节直接给出一套可复现的实验流程。我测试用的系统参数如下：H=5s，D=2，xd=1.81，x'd=0.3，T'd0=7.5s，系统等值电抗x'dΣ=0.5（包含变压器和线路），无穷大母线电压V∞=1.0。

5.1 数据生成：如何构造“真实”的动态轨迹

动态状态估计的仿真基准通常是用真值加噪声来模拟量测。我的做法是：先用ode45求解发电机三阶模型的连续方程，得到一条高精度的状态轨迹作为“真实值”ground truth，然后从这条轨迹里抽取离散时间样本，叠加高斯白噪声生成量测序列。

具体来说，仿真时长设为20s，扰动设置在t=5s，机械功率Pm从0.8 p.u.阶跃到1.0 p.u.。这个扰动会激起功角振荡，正好用来检验两种滤波器对动态过程的跟踪能力。采样周期取0.01s，每时刻给量测叠加1%的有功功率噪声和0.5%的电压幅值噪声，噪声协方差R就按这条设定。

滤波器的初始状态会略偏离真值（这是动态估计的常态，因为滤波启动时我们没有精确的初始状态），比如初始功角偏了0.05 rad、初始Δω偏了0.01 p.u.，这样能真实检验滤波器的收敛能力。

5.2 核心代码结构与关键函数

代码的整体结构采用脚本+函数句柄的方式：一个脚本负责生成真值、叠加噪声、调用滤波器、绘制结果；两个滤波器各写成一个函数文件；两个模型（状态函数f、量测函数h）写成匿名函数或子函数。下面是状态函数和量测函数的核心片段：

% 状态转移函数（离散化后的欧拉形式） f_state = @(x, u, dt, sys) [ x(1) + sys.wb * x(2) * dt; x(2) + ((sys.Pm - Pe(x, sys) - sys.D*x(2)) / (2*sys.H)) * dt; x(3) + ((sys.Efd - x(3) - (sys.xd - sys.xd1)*Id(x, sys)) / sys.Td0) * dt ];

量测函数Pe和Vt的计算建议直接写成嵌套函数，方便复用：

function pe = Pe_calc(x, sys) pe = (x(3) * sys.Vinf / sys.xds) * sin(x(1)); end

滤波器函数里，预测步和更新步严格分离，这样后续想换滤波器框架时只需要改函数入口和Sigma点生成逻辑，其他代码结构不用大改。

5.3 评估指标与结果分析

我用的评估指标是均方根误差RMSE和最大绝对误差MAE。RMSE看整体精度，MAE看最坏情况下的跟踪能力。跑完一次仿真，结果整理成类似下面这样的表：

从这个表里能看到几个很典型的结论：第一，UKF的RMSE大约是EKF的一半，尤其在扰动后的暂态段体现最明显；第二，单步耗时UKF确实更贵，但依然远低于采样周期，实时性完全没问题；第三，两者的稳态误差差别不大，差距全部发生在动态过程中。

画图时我习惯把真值、EKF估计值、UKF估计值叠在一张图里，再用一个子图单独画误差曲线。第一眼看上去三条线几乎重合，缩小到误差子图才能体现两者差异。这也是为什么很多论文的图看起来“EKF也挺准”——大部分场景下确实如此。

5.4 一个容易被忽略的问题：量测顺序对结果的影响

实际编程时，量测的接入顺序会影响立刻可见的观测效果。一个很常见的错误是在滤波初段就把所有量测并行接入，这会放大初始协方差确定不当的影响。我的做法是在初始阶段用较大的过程噪声Q让滤波器先自主收敛，等协方差矩阵收敛到合理水平后（通常几百个采样点之后）再将精确量测权重调大。这一招在实际调试中非常管用，能把滤波器从一个不太准的初值拉回来。

6. 调参踩坑实录：Q、R矩阵和其他看不见的细节

滤波算法的效果好坏，一半靠算法，一半靠参数。EKF和UKF对Q和R的敏感程度远超一般人的想象。

6.1 过程噪声Q矩阵：不是越大越好，也不是越小越好

Q矩阵在物理上代表我们对状态方程的信任程度。Q取小了，滤波器会过度信任模型，量测的修正作用被削弱，当模型本身有误差时容易发散；Q取大了，滤波器会跟着量测噪声乱跳，估计轨迹毛刺多，虽然“数据很新”，但精度反而下降。

我的经验是先用对角线定标：对功角任务，状态噪声量级按真值变化率的1%到5%来设。比如功角每步变化典型值0.1 rad，那Q(1,1)就取(0.001~0.005)²。这个经验值配合现场数据缩放，一般能拿到一个能用的起点。之后再针对具体场景做一次半自动调节，把量测仿真往真值上怼，观察到RMSE随Q变化的U形曲线，取谷底对应的值。

6.2 协方差矩阵的正定性维护

UKF对P矩阵的正定性要求比EKF苛刻得多，因为每次预测都要做Cholesky分解。我在一次仿真时发现，当功角接近π边界时，P矩阵的数值会突然半正定化，chol()直接报错。解决办法有两条：第一是给P加一个1e-12的Frobenius范数缩放对角阵，属于“保险丝”的加法；第二是改用sqrtm()函数求解平方根，它更稳定，代价是慢一点。工程上我建议两者都做：正常运行时用chol()，一旦抛异常就退化为sqrtm()，并把对小特征值截断到1e-10的量级。

6.3 滤波发散的监测与保护

动态估计中，滤波发散是最磨人、也最容易让结果完全没有可用性的问题。发散的典型症状是：新息序列（量测残差）不再服从零均值白噪声分布，而是出现系统性的偏差漂移。

我的调试习惯是在滤波器里加一个健康度指标，即每一步记录z - z_pred的Mahalanobis距离：d = (z - z_pred)' / (P_zz) * (z - z_pred)。正常情况下d的统计均值应该接近量测维数，如果连续几十步d都显著偏大，说明滤波器快不行了。遇到这种情况，常用的抢救措施是为Q加一个自适应调节因子：当新息连续偏大时主动涨Q，让滤波器重新“信任量测”拉回来。这个自适应策略在工程上很常见，但要注意别让它频繁触发，否则滤波器会变得过度活跃。

6.4 单位问题：p.u.系统和角度单位的统一

电力系统动态仿真里，角度量纲一会用弧度一会用度，这个细节看着不起眼，但能让你的结果差出百倍。我见过有人把δ的单位在量测方程里写成了度，而状态方程里用的是弧度，结果滤波器直接把量测当作完全错误的信息丢掉了。统一的做法是全程使用rad，仅在最后出图时转换为度。另外，转速偏差Δω的单位也要统一：如果是p.u.，需要乘上ωb才是rad/s。我的建议是状态空间里统一用rad/s，量测和绘图时再转换。

6.5 从单机系统向多机系统的扩展路径

这套单机实验的代码跑通后，向多机系统扩展的核心工作就是把状态方程和量测方程换成多机微分代数方程组，其余滤波框架完全不用动。多机系统的状态维数会涨到10到20维，此时EKF的雅可比矩阵解析推导会变得非常痛苦，数值差分也会因为维数上升而明显变慢；而UKF的Sigma点数量是2n+1，n=15时只需要31次函数求值，反而更有优势。我对那些想写多机动态估计的同学，一般都会建议直接用UKF作为主力滤波器，EKF用作对比基线就够了。

如果要进一步提高估计质量，可以考虑把UKF里的过程噪声Q做成随工况变化的自适应矩阵，或者引入交互多模型（IMM）来应对拓扑切换。这些方向都是动态状态估计领域比较活跃的研究点，从这篇基础实现往上扩展的路是通畅的。

回到Matlab工程本身，代码组织上我最后还想强调一点：不要把滤波算法和仿真数据生成混在同一个文件里。把系统参数、真值生成器、滤波器、评估指标分别封装成独立函数或脚本，后续换系统、换滤波器、换扰动场景时，你只需要改参数或换函数，不必重写全部代码。这也是我踩了两次“改一处牵全身”的坑之后总结出的最大教训。