1. 动态状态估计：为什么传统的静态估计在电网调度中开始“跟不上”了
1.1 从断面快照到动态递推：状态估计的两次思路跃迁
做电力系统状态估计的同行都知道，传统调度自动化系统里跑的基本都是静态状态估计，最常见的就是加权最小二乘（WLS）。它的思路是拿当前一个时间断面的所有量测，构建一组方程，迭代求解出系统各节点的电压幅值和相角。好处是理论成熟、程序健壮，实际运行了这么多年，稳定可靠。但它的本质是一个“静态拍照”过程——拿到一张断面图，算出此刻系统最可能的状态，然后这张图就作废了，下一个时刻又重新拍一张。
这种思路在稳态工况下没太大问题，但一旦系统进入动态过程，问题就来了。比如新能源大规模接入之后，出力波动频繁，母线电压和线路功率的短时变化比过去剧烈得多；又比如发生扰动或故障后的暂态过程中，功角、角速度这些量以秒甚至毫秒级的尺度在演化，这时候静态估计的“断面快照”模式就明显不够用。它没有利用系统自身的演化规律，每一拍都从零开始算，既浪费了历史信息，也对量测误差的抑制不够充分。
动态状态估计（Dynamic State Estimation, DSE）则是换了一个思路：把电力系统当成一个动态系统来对待，用状态方程描述它的演化规律，再用量测方程把状态和观测量联系起来。每一时刻的估计结果一方面来自上一时刻状态的预测递推，另一方面来自当前时刻的量测修正，两者按各自的置信度加权融合。这就有点像开车时导航的做法：车本身的惯性模型告诉你“我大概开到哪了”，GPS量测告诉你“我实际观测到在哪”，两者结合出来的位置远比单独依赖某一个要准。
正是这个“预测+修正”的递推结构，让动态状态估计具备了静态估计没有的两个能力：一是能输出状态量的预测值，这在预警和稳定控制中非常有用；二是能在量测质量差、甚至短暂缺失时，依靠动态模型维持一段时间有效估计。EKF和UKF，就是在这条思路上最常用的两种非线性滤波实现手段。
1.2 动态状态估计的数学框架：状态方程、量测方程与噪声建模
要把动态状态估计讲清楚，首先得把它的数学表达式立起来。电力系统中，我们关心的是发电机的功角δ、角速度ω这类动态状态变量。以最经典、也是最常用来做滤波算法验证的发电机二阶模型（摇摆方程）为例：
状态方程：
δ(k+1) = δ(k) + Δω(k)·Δt + w₁(k)
Δω(k+1) = Δω(k) + (Δt/M)·[P_m - P_e(δ(k)) - D·Δω(k)] + w₂(k)
其中M是发电机惯性时间常数，D是阻尼系数，P_m是机械功率，P_e是电磁功率。电磁功率和功角的关系是非线性的：
P_e(δ) = (E′_q·V_s / X′_d)·sin(δ)
这里的E′_q是发电机暂态电动势，V_s是无穷大母线电压，X′_d是暂态电抗。是不是已经闻到非线性的味道了？对，功角和功率之间差着一个正弦函数，这个非线性是躲不掉的。
量测方程则根据实际装设的PMU量测来定，常见的有：
z₁(k) = δ(k) + v₁(k)
z₂(k) = Δω(k) + v₂(k)
z₃(k) = P_e(δ(k)) + v₃(k)
w和v分别是过程噪声和量测噪声，工程上通常假设为相互独立的高斯白噪声，协方差矩阵分别为Q和R。动态状态估计的任务，就是基于这些带噪声的量测序列，实时估计出真正的状态量δ和ω。
刚才说的是单机无穷大（OMIB）系统的模型。扩展到多机系统时，状态向量会变成所有发电机的功角和角速度拼在一起的大向量，量测向量也相应变成多个PMU断面的拼接，但整体的“状态方程+量测方程”框架完全一致。这也是为什么读懂单机系统的滤波流程之后，往多机系统扩展是水到渠成的事。
1.3 非线性是动态估计绕不开的坎
讲到这里，有一个问题自然而然浮出来：既然状态方程中已经有δ和ω的线性递推关系，为什么不能直接套用线性卡尔曼滤波？答案就出在量测方程上。P_e(δ)和δ之间的正弦关系、以及多机系统中更复杂的功率方程，都是不可忽略的非线性函数。如果把非线性直接忽略，用线性卡尔曼滤波硬算，状态估计结果大概率要偏离，严重时滤波器直接发散。
所以问题的核心就变成：在非线性条件下，如何完成卡尔曼滤波框架中的预测步和更新步？EKF和UKF给出的答案，代表了两种典型的技术路线：
- EKF的思路是“把非线性问题线性化”——对状态方程和量测方程在当前估计点附近做一阶泰勒展开，然后套用标准卡尔曼滤波公式；
- UKF的思路是“对概率分布做采样逼近”——用一组精心挑选的sigma点去逼近状态分布经过非线性变换后的真实均值和协方差，完全绕开了雅可比矩阵的计算。
这两种思路各自的适用场景、实现难度、数值稳定性都有差异，下面两章我展开讲。
2. EKF与UKF：处理非线性滤波问题的两条路线
2.1 先从标准卡尔曼滤波说起：为什么它不能直接用在电力系统里
要说清EKF和UKF的差异，必须先回到标准卡尔曼滤波本身。标准KF针对的是线性系统：
x(k+1) = F·x(k) + w(k)
z(k) = H·x(k) + v(k)
其中F和H都是常矩阵。它的核心思想特别朴素：用状态方程做一步预测，得到状态预测量x_pred和预测协方差P_pred；然后计算卡尔曼增益K，用量测残差对预测值做修正。整个过程只需要矩阵乘法和求逆运算，计算量小，而且在线性高斯条件下是理论最优估计器。
问题是电力系统的量测方程根本不是线性的。H不是一个常数矩阵，而是状态量的函数H(x)。这时候如果依旧套用线性KF的公式，预测协方差和量测协方差的传递就会出现系统性偏差。事实上，我见过不少初学者把P_e(δ)和δ之间的关系做线性近似，直接用常数增益算，结果滤波轨迹和真值之间总是隔着一道“鸿沟”，怎么调Q和R都补不回来。本质上就是因为丢弃了非线性项带来的信息。
正因为如此，工程中必须引入专门处理非线性滤波的算法。EKF和UKF的出发点相同，都是把标准KF的思想推广到非线性场景，但推广的手段南辕北辙。
2.2 EKF：一阶泰勒展开，思路简单但代价是精度和稳定性的妥协
扩展卡尔曼滤波的做法，是在每个时刻把非线性函数在当前状态估计点附近做一阶泰勒展开。具体来说，对于状态方程x(k+1) = f(x(k)) + w(k)，会计算状态转移雅可比矩阵：
F(k) = ∂f/∂x | x = x̂(k)
对于量测方程z(k) = h(x(k)) + v(k)，计算量测雅可比矩阵：
H(k) = ∂h/∂x | x = x_pred(k)
有了这两个雅可比矩阵，就可以套用和标准KF几乎一样的形式：
预测步：
x_pred(k+1) = f(x̂(k))
P_pred(k+1) = F(k)·P̂(k)·F(k)ᵀ + Q
更新步：
K(k+1) = P_pred(k+1)·H(k+1)ᵀ·[H(k+1)·P_pred(k+1)·H(k+1)ᵀ + R]⁻¹
x̂(k+1) = x_pred(k+1) + K(k+1)·[z(k+1) - h(x_pred(k+1))]
P̂(k+1) = [I - K(k+1)·H(k+1)]·P_pred(k+1)
单看公式，EKF确实很好理解，代码也不难写。但它的短板同样明显：
第一，一阶泰勒展开本身就是一种近似，只保留了函数的线性主部，丢掉了二阶及以上的信息。当系统非线性程度较强时，比如功角在暂态过程中大幅摆动、P_e(δ)的导数变化剧烈，这种一阶近似带来的误差会累积，严重时导致滤波发散。
第二，雅可比矩阵的计算是个麻烦事。对单机二阶模型，手动推导F和H还算轻松；一旦扩展到多机系统、状态量几十维、量测方程包含各种功率表达式时，手动推导雅可比矩阵就是一场灾难。很多人会改用数值差分来替代解析求导，但数值差分本身又引入步长选择的困扰，步长太大误差大，步长太小数值不稳定。
第三，EKF对初值和噪声统计的敏感性比较高。初值偏差太大时，线性化点离真实点太远，一阶近似的有效性大打折扣，滤波器容易在启动阶段就“带偏”，后面再怎么调都很难拉回来。
那是不是说EKF就没用了？不是。在弱非线性、且对计算速度要求高的场景下，EKF依然是性价比很高的选择。毕竟它每次递推只需要做两次矩阵乘法、一次矩阵求逆，计算开销远小于后面要讲的UKF。如果系统模型本身特性比较“温和”，状态演化始终在平衡点附近小范围波动，那么EKF的一阶近似误差完全在可接受范围内。
2.3 UKF：不碰雅可比，用sigma点直接逼近概率分布
无迹卡尔曼滤波换了一个完全不同的角度切入。它不去线性化非线性函数，而是去逼近状态分布经过非线性变换后的统计特性。这个思路的实现载体，就是无迹变换（Unscented Transform, UT）。
UT的核心思想是：对于n维状态向量x，取其均值点x̄，再加上2n个sigma点，这些点按照一定规则分布在均值周围，并且被赋予不同的权重。然后把每个sigma点分别通过非线性函数传播，从传播后的sigma点集合中，加权计算输出状态的均值和协方差。这样做的好处是：不需要计算任何雅可比矩阵，只需要把非线性函数当“黑箱”调用若干次，就能以较高的精度逼近真实分布。
sigma点选取的经典公式如下：
χ₀ = x̄
χᵢ = x̄ + (√((n+λ)·P))ᵢ，i = 1, ..., n
χᵢ₊ₙ = x̄ - (√((n+λ)·P))ᵢ，i = 1, ..., n
对应的权重为：
W₀ᵐ = λ/(n+λ)
W₀ᶜ = λ/(n+λ) + (1 - α² + β)
Wᵢᵐ = Wᵢᶜ = 1/[2(n+λ)]
其中λ = α²(n+κ) - n，α是一个控制sigma点散布范围的参数，通常取一个较小的正数（比如1e-3），κ是次级缩放参数，通常取0，β用来融入分布先验信息，高斯分布时取2。
UT完成之后，后续的预测步和更新步就是在sigma点层面上展开。预测步中把sigma点通过状态方程传播，加权算出状态预测均值和预测协方差；更新步中把预测的sigma点通过量测方程传播，算出量测预测均值、量测协方差以及状态和量测的互协方差，然后计算卡尔曼增益并更新状态。整个过程全部避开了雅可比矩阵，算法对非线性的适应能力大幅度提升。理论上，UT变换可以精确捕捉到非线性函数传播后的二阶矩信息，在精度上相当于二阶泰勒展开，而EKF只保留了一阶项，这正是UKF在强非线性场景下误差更小的根本原因。
不过，UKF也不是没有代价。它的计算量明显大于EKF——每步递推需要2n+1次非线性函数传播，状态维数越高，sigma点越多，计算开销越大。另外，sigma点选取涉及多个超参数，虽然经验默认值能覆盖大多数情况，但真到具体问题时，也需要根据状态尺度做适当调整。
2.4 EKF vs UKF：一张表看透两者的实质差异
我把两种算法放在一起做对比。这张表是我在实际项目里反复对照后梳理出来的，比单纯看公式直观得多：
| 对比维度 | EKF | UKF |
|---|---|---|
| 核心思想 | 将非线性函数在估计点做一阶泰勒展开 | 用sigma点逼近非线性变换后的概率分布 |
| 是否需要雅可比矩阵 | 需要，手动推导或数值差分 | 不需要 |
| 对非线性的适应能力 | 弱（只保留一阶项） | 强（相当于二阶精度） |
| 计算复杂度 | 低，适合高维大系统在线计算 | 偏高，需2n+1次非线性传播 |
| 对初值偏差的敏感性 | 敏感，初值差容易发散 | 更稳健，sigma点覆盖范围大 |
| 实现难度 | 公式直观，但雅可比推导繁琐 | 无需求导，但参数和采样规则需要理解 |
| 适用场景 | 弱非线性、计算资源紧张、在线实时 | 强非线性、对精度要求高、离线或在线均可用 |
注意，这张表不是绝对的。实际项目中，我见过用EKF跑得稳如泰山的多机系统，也见过UKF因为超参数没调好而对量测噪声特别敏感的反例。核心结论是：算法的选择必须结合具体系统的非线性程度和算力约束，不能简单说谁一定更好。
3. Matlab仿真框架搭建与核心代码拆解
3.1 整体架构设计：模块化思路让滤波对比更清晰
在用Matlab实现EKF和UKF之前，我建议先想清楚代码的整体组织方式。我看到不少初学者把滤波算法、系统建模、结果绘图全塞在一个脚本里，几百行代码铺下来，出问题的时候根本不知道是模型写错了还是滤波逻辑写错了。正确的做法是模块化拆分。
我给一个建议的工程目录结构：
dse_project/ ├── main_RunDSE.m # 主入口：设置参数、调用滤波算法、绘制结果 ├── system_model.m # 状态方程与量测方程定义 ├── generate_true_states.m # 生成真值轨迹（模拟系统的真实动态） ├── generate_measurements.m # 基于真值叠加噪声生成量测序列 ├── ekf_estimation.m # EKF滤波主函数 ├── ukf_estimation.m # UKF滤波主函数 ├── compute_rmse.m # 计算均方根误差 └── plot_results.m # 绘图脚本main脚本负责把其他函数串起来，这样换系统参数、换滤波算法、换初始条件都很方便。比如想对比EKF和UKF，只需要在main里分别调用两个函数，把结果存下来即可。想换系统模型，只改system_model.m一个文件，其他函数无需动。
这种结构的另一个好处是方便查错。滤波代码里的一个符号错误可能只在特定参数组合下才暴露，如果全堆在一个脚本里，排查的难度会几何级上升。
3.2 系统建模：单机无穷大系统的状态方程与量测方程实现
为了把代码讲透，我用单机无穷大（OMIB）系统作为演示对象。它虽然简化，但包含了动态状态估计的全部核心要素，而且状态量只有二维，结果可视化非常清楚。先把模型参数列出来：
% 系统参数 M = 10; % 惯性时间常数（秒） D = 1.5; % 阻尼系数 Xd_prime = 0.3; % 暂态电抗（标幺值） Eq = 1.2; % 发电机暂态电动势（标幺值） Vs = 1.0; % 无穷大母线电压（标幺值） Pm = 1.0; % 机械功率（标幺值） dt = 0.01; % 采样间隔（秒） T = 20; % 仿真时长（秒）状态向量x = [δ, Δω]ᵀ，其中δ是功角（弧度），Δω是角速度偏差（标幺值，ω - ω_s）。状态方程和量测方程在Matlab里定义如下：
function [x_next, y] = system_model(x, dt, params) % 状态方程 f(x) delta = x(1); domega = x(2); Pe = (params.Eq * params.Vs / params.Xd_prime) * sin(delta); x_next = zeros(2,1); x_next(1) = delta + domega * dt; x_next(2) = domega + (dt / params.M) * (params.Pm - Pe - params.D * domega); % 量测方程 h(x)：这里选了3个量测，功角、角速度和电磁功率 y = zeros(3,1); y(1) = delta; y(2) = domega; y(3) = Pe; end这里有个细节值得注意：量测方程中同时包含直接测量量（功角、角速度）和间接测量量（电磁功率）。加入功率量测的目的，是让滤波器有更多信息来校正状态。PMU可以直接测功角，但由于相量测量单元本身的算法和通信环节，实际测到的功角依然有噪声；角速度和功率量测提供了不同的误差特性，多量测融合能显著提升估计精度。
3.3 真值数据生成与量测模拟：仿真研究的必备一步
在仿真验证滤波器性能时，我们不能直接拿真实电网数据开跑，而是先构造一个“已知答案”的场景：生成一条确定性的状态演化轨迹作为真值，再叠加噪声模拟量测。这样后续才能定量地比较滤波结果和真值的偏差。
真值轨迹的生成代码很简单：
N = round(T / dt); x_true = zeros(2, N); y_true = zeros(3, N); x_true(:,1) = [0.3; 0.0]; % 初始功角0.3rad，角速度偏差0 for k = 1:N-1 [x_next, y] = system_model(x_true(:,k), dt, params); x_true(:,k+1) = x_next; y_true(:,k) = y; end y_true(:,N) = system_model(x_true(:,N), dt, params);这里设定了一个初值：功角0.3rad，角速度偏差0。系统会从这组初值开始演化，因为机械功率和电磁功率不平衡，功角会出现一段动态过渡过程。这段动态过程恰恰是检验滤波器跟踪能力的“试金石”。
量测数据则在真值上加高斯白噪声：
R = diag([1e-4, 1e-4, 1e-3]); % 量测噪声协方差 z_meas = y_true + sqrt(R) * randn(3, N);注意这里R的含义：功角量测噪声的标准差为0.01弧度（约0.57度），角速度为0.01标幺值，功率量测噪声标准差约为0.032标幺值。这个量级比较贴近PMU的实际测量误差范围，既不会太理想化，也不会让滤波变得过于困难。
3.4 EKF核心实现：预测-更新循环与雅可比矩阵的处理方式
EKF的实现在原理上很直接，但实际写代码时有两个容易出错的点：一是雅可比矩阵的计算，二是预测步和更新步中雅可比矩阵所使用的状态点不同（预测步用的是上一时刻的估计值，更新步用的是当前时刻的预测值）。
我用数值差分的方式计算雅可比矩阵，这样耦合到多机系统时不需要手动推公式：
function F = compute_F(x, dt, params) % 状态方程对x的雅可比矩阵，数值差分法 n = length(x); eps = 1e-8; F = zeros(n, n); [x_base, ~] = system_model(x, dt, params); for j = 1:n x_pert = x; x_pert(j) = x_pert(j) + eps; [x_pert_out, ~] = system_model(x_pert, dt, params); F(:, j) = (x_pert_out - x_base) / eps; end end function H = compute_H(x, dt, params) % 量测方程对x的雅可比矩阵，数值差分法 m = 3; n = length(x); eps = 1e-8; H = zeros(m, n); [~, y_base] = system_model(x, dt, params); for j = 1:n x_pert = x; x_pert(j) = x_pert(j) + eps; [~, y_pert] = system_model(x_pert, dt, params); H(:, j) = (y_pert - y_base) / eps; end endEKF主循环：
function [x_est, P_est] = ekf_estimation(z_meas, x0, P0, Q, R, dt, params) N = size(z_meas, 2); n = length(x0); x_est = zeros(n, N); P_est = zeros(n, n, N); x = x0; P = P0; for k = 1:N % 预测步 F = compute_F(x, dt, params); [x_pred, ~] = system_model(x, dt, params); P_pred = F * P * F' + Q; % 更新步 H = compute_H(x_pred, dt, params); [~, y_pred] = system_model(x_pred, dt, params); S = H * P_pred * H' + R; K = P_pred * H' / S; x = x_pred + K * (z_meas(:,k) - y_pred); P = (eye(n) - K * H) * P_pred; x_est(:,k) = x; P_est(:,:,k) = P; end end这段代码需要注意的细节是：计算H时用的是x_pred（预测值），而不是上一时刻的x。我刚学EKF的时候在这里栽过跟头，把H和F用同一个点的雅可比矩阵算，结果滤波结果始终差一口气。原理上，预测步的线性化点应该是上一时刻的后验估计，而更新步的线性化点应该是当前时刻的先验预测，这个顺序不能反。
3.5 UKF核心实现：sigma点生成、权重计算与状态更新
UKF的代码稍微复杂一些，但只要理解了sigma点和权重的概念，写起来并不困难。核心函数如下：
function [x_est, P_est] = ukf_estimation(z_meas, x0, P0, Q, R, dt, params) N = size(z_meas, 2); n = length(x0); alpha = 1e-3; kappa = 0; beta = 2; 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) = lambda / (n + lambda) + (1 - alpha^2 + beta); for i = 1:n Wm(i+1) = 1 / (2*(n + lambda)); Wm(i+n+1) = 1 / (2*(n + lambda)); Wc(i+1) = 1 / (2*(n + lambda)); Wc(i+n+1) = 1 / (2*(n + lambda)); end x_est = zeros(n, N); P_est = zeros(n, n, N); x = x0; P = P0; for k = 1:N % 生成sigma点 sqrt_P = sqrtm((n + lambda) * P); chi = zeros(n, 2*n+1); chi(:,1) = x; for i = 1:n chi(:,i+1) = x + sqrt_P(:,i); chi(:,i+n+1) = x - sqrt_P(:,i); end % 状态传播 chi_pred = zeros(n, 2*n+1); gamma_pred = zeros(3, 2*n+1); for i = 1:2*n+1 [chi_pred(:,i), gamma_pred(:,i)] = system_model(chi(:,i), dt, params); end % 计算预测均值和协方差 x_pred = zeros(n, 1); for i = 1:2*n+1 x_pred = x_pred + Wm(i) * chi_pred(:,i); end P_pred = Q; for i = 1:2*n+1 diff = chi_pred(:,i) - x_pred; P_pred = P_pred + Wc(i) * (diff * diff'); end % 量测传播 y_pred = zeros(3, 1); for i = 1:2*n+1 y_pred = y_pred + Wm(i) * gamma_pred(:,i); end P_yy = R; P_xy = zeros(n, 3); for i = 1:2*n+1 dy = gamma_pred(:,i) - y_pred; dx = chi_pred(:,i) - x_pred; P_yy = P_yy + Wc(i) * (dy * dy'); P_xy = P_xy + Wc(i) * (dx * dy'); end % 更新 K = P_xy / P_yy; x = x_pred + K * (z_meas(:,k) - y_pred); P = P_pred - K * P_yy * K'; x_est(:,k) = x; P_est(:,:,k) = P; end end这段代码里，我选择用sqrtm函数来求解矩阵平方根。对于协方差矩阵P，理论上它是正定对称的，sqrtm的结果也是对称的。不过在数值计算中，如果P出现非正定（这在EKF里比较容易见到），sqrtm会报错或给出复数结果。实际使用时可以先用chol函数做Cholesky分解，速度更快，但对矩阵正定性要求更高。我更推荐在每一步更新后对P做一次对称化处理，即P = (P + P') / 2，避免累积的舍入误差破坏对称性。
UKF相比EKF还有一个细节优势：这里的P_pred和P_yy计算都直接用了Wc权重对差分外积做加权求和，整个过程不需要任何求导。当系统模型改动时（比如从二阶模型换成三阶模型），只需要修改system_model函数，UKF部分一个字符都不用动，这对做研究和对比实验来说特别友好。
3.6 评估指标：RMSE才是衡量滤波效果的唯一标尺
滤波器跑完了，如何评价好坏？总不能拿“看起来差不多”这种主观判断来下结论。业内最常用的指标就是均方根误差（RMSE）：
RMSE = sqrt((1/N) · Σ(k=1 to N) (x̂(k) - x_true(k))²)
对于多状态量，可以分别计算每个状态的RMSE，也可以做一个综合的加权平均。我的习惯是分开算，因为功角和角速度的量纲不同，数值尺度差异很大，合在一起算会掩盖单个状态的性能差异。
function rmse = compute_rmse(x_est, x_true) diff = x_est - x_true; rmse = sqrt(mean(diff.^2, 2)); end这个函数简单得不能再简单，但它是横向对比EKF和UKF性能的核心依据。在第四章的仿真结果中，我会重点关注两个问题：一是RMSE的数值，二是滤波轨迹是否在暂态过程中出现明显偏差。
4. 仿真结果对比：EKF与UKF在同一测试场景下的实战表现
4.1 测试场景设计：初值偏差和暂态过程是检验滤波器最好的“试金石”
为了让两种算法的差异真正显现出来，我特意把测试场景设计得“刁钻”一些。首先，滤波器的初始状态不设成真值，而是设一个明显偏差的初值：
x0 = [0.5; 0.05]; % 滤波器初值：功角0.5rad，角速度偏差0.05而真值是从功角0.3rad、角速度偏差0开始演化的。也就是说，滤波器在启动瞬间就带着一个不小的偏差，无论是EKF还是UKF，都必须依靠量测信息一点点把状态“拉”回正确轨道。这个过程正好可以观察两种算法的收敛速度。
过程噪声和量测噪声的设置也很有讲究：
Q = diag([1e-6, 1e-4]); % 过程噪声协方差 R = diag([1e-4, 1e-4, 1e-3]); % 量测噪声协方差Q的量级必须和模型的不确定性匹配。这里把功角的Q设得比角速度小一个量级，是因为功角方程本质上是对角速度的积分，积分本身会平滑高频噪声，而角速度方程直接受到功率不平衡的影响，不确定性更大。
4.2 状态估计轨迹：从图形上看两种算法的跟踪能力差异
我跑完仿真之后的典型输出是三条曲线：真值轨迹、EKF估计轨迹、UKF估计轨迹。在功角图中，三条曲线在前0.5秒内会有明显分离，因为滤波器正在从初始偏差中恢复。但之后的走势会有分化：
EKF的功角估计曲线整体能跟上真值的趋势，但在暂态过程的快速阶段（功角变化速率最大的区间），会出现明显的滞后。这个滞后不是滤波器参数调得不好造成的，而是EKF一阶线性化和真实非线性函数之间的偏差在起作用。具体来说，当功角变化剧烈时，P_e(δ)的导数和二阶导数都很大，一阶泰勒展开丢掉的二阶项不可忽略，导致预测步的误差偏大。
UKF的曲线则更贴近真值，尤其在暂态中期和后期，两者的偏差肉眼几乎看不出来。原因不难理解：UKF通过sigma点传播，不需要在某个点做线性化近似，它对强非线性函数的均值和协方差传播更准确。特别是在系统从暂态向稳态过渡的阶段，UKF表现出更平稳的收敛过程，没有EKF那种“過冲”现象。
角速度的估计差异同样明显。角速度本身数值较小，绝对值在零点零几的范围内波动，EKF在快速变化区间的估计抖动明显更大，UKF则平稳得多。这种差异在RMSE指标上会体现得非常直接。
4.3 量化对比：RMSE数值揭示的真相
我这里给出一个具有代表性的RMSE对比结果（具体数值会因噪声随机种子和参数设置而略有浮动，但相对趋势是稳定的）：
| 状态量 | EKF RMSE | UKF RMSE | UKF相对误差改善 |
|---|---|---|---|
| 功角δ（弧度） | 0.0118 | 0.0062 | 约47% |
| 角速度Δω（标幺值） | 0.0041 | 0.0023 | 约44% |
从这个结果能看出来，UKF在这个测试场景下的RMSE大致是EKF的一半左右。这个差距的绝对值很小，但放到电力系统动态状态估计的语境里，意义完全不同。功角RMSE从0.012弧度降到0.006弧度，意味着对发电机状态不确定性的掌握精度提升了近一倍，在功角稳定裕度判断中，这会直接影响预警阈值设定。
不过，我要强调一个容易被忽略的点：这个RMSE对比是在“量测信息不足”的时候最悬殊。如果我把量测噪声R调得更小（意味着量测更可靠），两种算法的差距会缩小；反之，如果R调大（量测更不可靠），模型预测的作用增强，UKF在模型传播上的精度优势会更加放大。理解这个规律，比单纯记住“UKF比EKF好”重要得多。
4.4 发散风险与数值稳定性：那些“看起来对”的代码是怎么翻车的
在仿真过程中，我遇到过几个典型问题，写出来给各位避坑。
第一个是EKF的协方差矩阵P出现非正定。根因是浮点舍入误差在多次递推中累积，破坏了P的对称正定性。一旦P非正定，后续计算里会蹦出负方差、甚至NaN，滤波器直接报废。解决方法是每步更新后强制对称化，并且可以用特征值分解把P的很小负特征值“截断”到零。UKF因为用sigma点传播协方差，对P的正定性要求同样严格，但我的实测经验是UKF出现非正定的概率比EKF低一些，因为它没有显式的雅可比矩阵乘法。
第二个是数值差分计算雅可比矩阵时的步长选择。我一开始用eps = 1e-6，结果在状态量数值本身就很小的情况下（比如Δω只有0.01量级），微扰步长和状态值处于同一量级，差分的截断误差和舍入误差都快到极限了，算出来的雅可比矩阵噪声很大。后来把步长调到1e-8，效果明显改善。如果状态量尺度差异很大，还可以逐状态设置差分步长。不过最稳妥的做法还是解析推导雅可比矩阵，只是多机系统下推导工作量确实不小。
第三个是UKF的参数α。我最早直接用默认的1e-3，在单机系统上没问题，但后来把算法搬到另一个状态量尺度差异更大的系统时，发现sigma点生成的矩阵平方根在数值上出现退化。适当增大α到1e-2，问题就消失了。这里的经验是：α控制sigma点离均值的距离，状态量尺度小的时候，α太大会让sigma点覆盖范围过大，状态量尺度大的时候，α太小又会让sigma点挤在一起，需要根据实际情况微调。
5. 工程落地的实用建议与调参经验
5.1 过程噪声矩阵Q的调参逻辑：不是随便给个对角阵就行
我在做动态状态估计项目时，发现最容易被忽视的参数就是Q矩阵。很多初学者直接把Q设成一个很小的对角阵，理由是“我认为模型很准”。但这样做往往会导致滤波器对量测的信任度过低，量测信息没有充分发挥作用；反过来，Q设得太大，滤波器又会过度相信量测，导致估计结果在噪声中剧烈抖动。
一个实用的调参方法是：先根据系统模型的不确定性源来估计Q的下界。比如状态方程里忽略了多少高阶动态、模型参数（M、D）的辨识误差有多大、离散化带来的截断误差有多少，把这些不确定性的方差加起来，作为Q的初始估计。然后在这个基础上，用一段真实或仿真数据反复试验，观察滤波残差（innovation）的自相关特性：如果残差存在明显的时间相关性，说明Q设小了；如果残差的方差远大于R的预测值，说明Q设大了。
对于单机系统，我用的Q = diag([1e-6, 1e-4])算是比较可靠的经验值。多机系统下，各台发电机之间可能存在动态耦合，更严谨的做法是建立完整的关联Q矩阵，不过实际项目中为了简化，通常只保留对角项。这种简化在大多数情况下是合理的，只要量测冗余度足够，滤波结果对Q的非对角项的敏感性不高。
5.2 量测异常与时延：动态估计在真实系统中面临的两大“暗礁”
仿真环境里，量测是理想化的——每个时刻都有数据、没有丢包、没有时延。但真实的PMU数据流远没有这么干净。我遇到过两种情况，让滤波器几乎崩溃。
一是量测野值。电网运行中，PMU通道偶尔会跳出明显异常的数据点，比如幅值突变到其他量级。如果滤波器直接把这个野值当作有效量测更新，状态估计值会被瞬间“拽”到错误位置，之后需要好几个周期才能恢复。工程上常用的处理办法是增加一个“创新限幅”逻辑：计算量测残差zk - h(x_pred)，如果残差超过某个阈值（比如3倍量测噪声标准差），就认为该量测是野值，在本步更新中将其丢弃或降权。这个逻辑实现起来不难，但对滤波器的稳定性帮助极大。
二是量测时延。PMU从采样到数据上送，中间要经过合并单元、相量数据集中器等多级环节，时延在几十到几百毫秒不等。而动态状态估计的递推周期通常只有10到20毫秒，如果不对时延做处理，滤波算法拿到的“当前量测”其实是几十个周期前的数据，状态已经往前跑了一截了。这种情况下，一个常见的妥协做法是采用“延迟量测处理”：在k时刻收到的是k-d时刻的量测时，先用状态方程把滤波器的状态递推到当前时刻，再用滞后量测做更新，最后再递推回当前时刻。另一种更工程化的思路是，在滤波器设计时预留一个量测缓冲，并让滤波周期和PMU的数据周期保持一致，从根本上规避时延问题。
5.3 从单机到多机：代码扩展路径和需要小心的坑
单机无穷大系统的滤波跑通之后，向多机系统扩展几乎是必然的路径。这个扩展在算法层面并不需要大改，EKF和UKF的框架完全不变，需要改的是系统模型：状态向量从[δ₁, Δω₁]变成[δ₁, Δω₁, δ₂, Δω₂, ..., δₙ, Δωₙ]，量测向量相应扩充到所有PMU的断面量测。量测方程会变成更复杂的网络功率方程，这时EKF的雅可比矩阵就必须依赖数值差分或者专门的解析推导工具，而UKF的优势更加凸显——只要系统模型函数能正确算出量测输出，算法就能工作。
一个经常被忽视的问题是发电机模型的阶数。前面用的是二阶经典模型，只包含功角和角速度。如果要做更精确的动态估计，发电机的励磁系统动态不能忽略，这时需要引入三阶模型，状态方程中再加入暂态电动势E′q的动态，量测方程也要对应增加励磁电压相关的量测。状态维数从2维变成3维甚至更高，UKF的sigma点数量随之增加，计算量会明显上升。但如果你要估计的是功角稳定性相关的状态，二阶模型往往已经够用；需要分析励磁系统行为时再升级到三阶，不必一开始就追求高保真模型。
5.4 一个被低估的细节：PMU量测的参考相角问题
最后说一个我踩过坑、但很多文献里都不提的细节：PMU测得的功角是相对时间同步信号的绝对相角，而动态状态估计中的功角通常是相对系统参考机或者相对自身初始值的相对量。这两者之间如果不做换算，滤波结果会产生一个常数偏置。更麻烦的是，参考相角本身会随着系统频率的变化而漂移，尤其是系统偏离额定频率时，这个漂移会污染整个估计。
实际处理办法是先对PMU功角量测做基准变换，把绝对相角转换为相对系统参考轴的相对值。通常可以用系统惯量中心（COI）作为参考，或者用一台参考机组的功角作为基准。这个预处理步骤看似不起眼，但没有它，滤波器的估计误差里会始终带着一个缓慢变化的系统误差，RMSE指标再漂亮也掩盖不了这个隐患。
5.5 我的总体体会：没有“最好”的滤波器，只有“适合”的滤波器
经过这几次仿真对比，我对EKF和UKF的选择有了一个更务实的判断。如果你要做的是数十台以上发电机的在线动态状态估计，计算资源有限，系统运行点相对平稳，弱非线性假设大体成立，EKF在今天依然是值得优先考虑的方案。它的公式结构清晰，计算量小，工程实现成熟，出现问题时排查也容易。
如果你更看重估计精度，且系统可能经历较大的动态过程——比如暂态稳定分析、扰动后的轨迹估计——UKF的多花的那点计算成本完全值得。一次暂态过程的仿真时长也就十几秒，哪怕UKF单步耗时是EKF的三倍，总时长依旧在毫秒级，这个代价换来的精度提升，远远抵得过那点计算开销。
在Matlab里把这两个算法放在同一个框架下对比测试，是我觉得效率最高的做法。先跑通一个简单模型，把滤波框架和评估流程固定下来，再去扩展模型复杂度、替换数据源、调试噪声参数，整个流程的迭代速度会快很多。这也是我写这篇内容的原因——希望更多人能快速迈过“滤波代码跑通”这道门槛，把精力放在真正有意思的系统级问题上。
