最近在做电力系统动态状态估计相关的仿真项目，把扩展卡尔曼滤波（EKF）和无迹卡尔曼滤波（UKF）在Matlab里的实现从头到尾完整捋了一遍。这套东西我在单机无穷大系统、多机系统上都跑过，也和实验室几个师弟交流过踩坑经历，发现很多人并不是卡在算法原理上，而是卡在“代码怎么组织、参数怎么设、跑出来怎么判断好坏”这些实战细节上。所以这篇文章就把整套实现方案、核心代码、调参过程和排错经验都写清楚，想用Matlab做电力系统动态状态估计的同学可以直接参考。

这篇内容适合正在做电力系统动态状态估计论文的研究生、想用PMU量测数据做在线状态估计的工程师，也适合第一次接触非线性卡尔曼滤波、想在Matlab里把EKF和UKF跑通并对比性能的初学者。你不需要一开始就弄懂所有数学推导，跟着代码把流程跑起来，再回头看原理，理解成本会低很多。

1. 动态状态估计要解决什么问题，为什么选EKF和UKF

1.1 从静态估计到动态估计：多了一个时间维

传统电力系统状态估计，业内叫静态状态估计，本质上是一个加权最小二乘问题。它把某个时间断面上所有SCADA量测拿过来，以节点电压幅值和相角为状态量，离线迭代求解一个最优断面。好处是理论成熟、实现简单，坏处也很明显：它把每一个时间断面当成孤立的，上一时刻怎么变化、系统的动态方程是什么样的，统统不管。

带来的实际问题是：量测冗余度不够时，断面估计可能完全失真；坏数据检测需要额外花功夫；而且SCADA量测刷新率低，秒级甚至分钟级，根本追不上系统动态。真正需要在线监视发电机功角、转速等动态状态的时候，静态估计就力不从心了。我曾经在一个含风电接入的算例里对比过，故障后半个周波内，静态估计给出的功角可能还停在故障前的位置，这对广域保护决策来说基本没用。

动态状态估计就不一样，它引入了系统的状态转移方程，把状态估计变成一个递推问题。上一时刻的状态通过发电机运动方程、电磁暂态方程预测到当前时刻，再用当前量测去修正。这样做的收益是：量测瞬时缺失时，预测值仍然能给出一段可用的状态轨迹；量测有坏数据时，滤波器天然有平滑作用；更重要的是，它能给出状态的真实动态变化过程，这直接服务于动态安全评估、广域控制、参数辨识、以及新能源并网后的惯量监测。

动态状态估计的状态量一般取发电机的功角、角速度偏差、暂态电动势等物理量。这些量没法直接用表计测量，但可以通过PMU测得的功角、电压幅值等间接推算出来。

1.2 非线性滤波选型：EKF、UKF之外的取舍

状态方程确定以后，下一个问题是：用什么滤波器来递推。线性卡尔曼滤波在电力系统里几乎没有直接用武之地，因为发电机动态方程是强非线性的，比如电磁功率Pe与功角δ之间的关系是正弦函数，摆动方程里这个非线性躲都躲不掉。所以必须在非线性滤波框架里选型。

EKF的思路最直白：在每一个时刻把非线性函数做一阶泰勒展开，只保留线性项，然后套用标准卡尔曼滤波的预测-更新流程。它需要显式计算雅可比矩阵，实现简单、计算量小，是很多人入门非线性滤波的第一站。但它的短板也出在线性化上：系统非线性越强、仿真步长越大，一阶近似误差越大；遇到强扰动或者状态突变，滤波结果可能出现明显滞后甚至发散。

UKF走的是另一条路线，不线性化函数本身，而是用一个精心设计的Sigma点集合去“采样”状态分布，让每个点通过非线性函数，再从传回来的点集里重新统计均值和协方差。它在数学上能捕捉到三阶矩的精度，不需要求雅可比矩阵，对强非线性系统的适应能力明显好于EKF。

为什么不直接上粒子滤波？粒子滤波处理非线性和非高斯问题的能力最强，但代价是需要大量粒子，计算量比UKF高两三个数量级，对于需要在线实时递推的动态状态估计来说，工程上很难接受。我在实际对比中也验证了这一点：EKF调得好时稳态跟踪没什么大问题，但故障暂态期间UKF明显收敛更快、误差更小；而粒子滤波在同样场景下虽然精度最高，仿真一个10秒算例要跑几分钟，在线应用根本不现实。所以如果你的研究目标不是专门做粒子滤波算法改进，EKF和UKF是最优的折中选择。

2. 算法核心拆解：EKF和UKF的数学逻辑与适用边界

2.1 EKF：非线性函数的“线性化近似”

EKF处理的系统可以写成：

x(k+1) = f(x(k), u(k)) + w(k) z(k) = h(x(k)) + v(k)

其中w为过程噪声，v为量测噪声，分别假设为高斯白噪声，协方差矩阵为Q和R。

预测步做两件事。第一是用状态方程直接积分，得到先验状态估计；第二是把状态方程在上一时刻估计值处做一阶泰勒展开，计算雅可比矩阵F，再通过F把上一时刻的协方差P变换到预测协方差。公式是：

x_pred = f(x_hat(k)) P_pred = F * P * F' + Q

更新步也要算量测方程的雅可比矩阵H，然后计算卡尔曼增益K：

K = P_pred * H' * (H * P_pred * H' + R)^(-1) x_hat(k+1) = x_pred + K * (z(k+1) - h(x_pred)) P(k+1) = (I - K*H) * P_pred

这里的关键问题有两个。第一个是雅可比矩阵F和H必须每一步重新算，因为工作点一直在变。第二个是雅可比矩阵的推导很容易出错，状态方程里只要有一项偏导写错，滤波器轻则精度下降，重则直接发散。我去实验室看过不少同学的代码，EKF结果诡异的基本都是这个原因。

一个非常实用的替代方案是：先用数值差分替代解析雅可比，把主流程跑通，确认滤波逻辑没有问题，再花时间推导解析式来加速。这个思路我强烈建议新手采用。数值雅可比本质上就是用中心差分近似偏导数，对状态量每个分量给一个小扰动，计算函数输出的变化率。只要扰动步长取值合适，精度完全够用。

2.2 UKF：Sigma点采样下的统计逼近

UKF的核心是无迹变换（Unscented Transformation）。它不直接对非线性函数做近似，而是对状态的分布做近似。假设状态向量是n维的，均值为x，协方差为P，那么生成2n+1个Sigma点：

X^(0) = x X^(i) = x + (sqrt((n+lambda)*P))_i, i = 1,...,n X^(i+n) = x - (sqrt((n+lambda)*P))_i, i = 1,...,n

其中sqrt((n+lambda)*P)表示对协方差矩阵做Cholesky分解或者特征值分解后取矩阵平方根的第i列。lambda是缩放参数，由alpha、beta、kappa三个参数决定：

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

alpha控制Sigma点围绕均值的散布程度，通常取1e-3到1之间的值；beta用来融进状态分布的先验信息，高斯分布下取2最优；kappa是次级缩放参数，一般取0或者3-n。权重按下式计算：

Wm^(0) = lambda / (n + lambda) Wc^(0) = lambda / (n + lambda) + (1 - alpha^2 + beta) Wm^(i) = Wc^(i) = 1 / (2*(n + lambda)), i = 1,...,2n

得到Sigma点以后，每个点都通过状态方程传递到下一时刻，再按权重重新合成预测均值和预测协方差。量测更新也做同样的操作：所有Sigma点通过量测方程，得到量测预测的均值和协方差，再算状态与量测的交叉协方差，最后计算增益并更新状态。整个过程不需要求任何雅可比矩阵。

我用一个类比来理解EKF和UKF的区别：EKF就像你拿到一张弯曲的地图，非要在当前所在位置用一把直尺把地图抹平，然后用直线去量距离；UKF则是派出一小队侦察兵往不同方向走，每个人回来后汇报海拔和坐标，再把这些信息加权汇总。显而易见，地形越复杂，侦察兵方案越可靠。

UKF的数值实现有一个绕不开的坑：每一步生成Sigma点都需要对协方差矩阵做平方根分解，而协方差矩阵在递推过程中可能因为舍入误差失去正定性，导致chol函数报错。我的处理办法是给协方差矩阵加一个很小的对角修正项，比如1e-8乘以单位阵，保证分解能正常进行。这个细节很多教程不会写，但跑起来迟早会遇到。

2.3 两者对比：选型逻辑和工程判断

实际工程里怎么选？我的建议是：如果你的状态方程非线性不严重，比如主要靠量测方程修正，那EKF足够，没必要为一个模型精度问题引入更多计算量；但如果你的算例里扰动大、状态变化剧烈，或者你根本推不出解析雅可比矩阵，直接上UKF，省心且效果好。做研究写论文的时候，两个都实现并做对比是最稳妥的，审稿人也喜欢看到这种比较。

3. Matlab实现：从系统建模到主循环代码

3.1 三阶发电机模型与量测方程设计

为了让代码既能说明问题又不过度复杂，我选了单机无穷大系统（SMIB）配三阶发电机模型。状态向量取三个量：

x(1) = delta 功角（rad） x(2) = dw 角速度偏差（标幺值，实际工程中也可用rad/s） x(3) = Eq 暂态电动势（标幺值）

状态方程写成：

d(delta)/dt = wb * dw d(dw)/dt = (Pm - Pe - D*dw) / M d(Eq)/dt = (Efd - Eq - (xd - xdp)*Id) / Tdo

其中wb = 2pif0是同步角速度基准，M = 2H是惯性时间常数，D是阻尼系数，Pm是机械功率，Pe是电磁功率，Efd是励磁电动势，Tdo是励磁绕组时间常数。电磁功率和电流按下式计算：

Id = (Eq - Vinf*cos(delta)) / xSum Iq = Vinf*sin(delta) / xSum Pe = Eq*Vinf*sin(delta) / xSum

这里的xSum是暂态电抗与线路电抗之和。量测方程取PMU能提供的两个量：功角delta和机端电压幅值Vt，其中Vt由d轴和q轴分量合成：

Vd = xq * Iq Vq = Eq - xdp * Id Vt = sqrt(Vd^2 + Vq^2)

仿真参数我建议这样取：

初始运行点通过稳态条件反推：先给定Eq0，再根据Pm与Pe相等的条件求出delta0。仿真步长取0.01s，对应100Hz的PMU量测频率，这个步长下Euler积分已经有点吃力，建议后续换成RK4。

3.2 数据生成与初始化：先有真值，再有量测

动态状态估计的仿真逻辑和真实系统是反过来的：在仿真环境里，我们先用高精度数值积分（比如ode45）跑出一条“真实轨迹”，再往真实轨迹上叠加高斯噪声模拟PMU量测，最后拿含噪声的量测去喂给EKF和UKF，并把它们的估计结果与真实轨迹做对比。这个流程很多人都知道，但实际操作中有一个细节常被忽略：千万不要把无损轨迹直接当量测输入滤波器，那样你会得到一个低得离谱的估计误差，而且完全检验不出算法的抗噪声能力。

量测噪声的生成也有一点讲究。角度量测噪声要换算成弧度，PMU的角度误差通常在0.01度到0.1度之间，对应弧度范围大约是1.7e-4到1.7e-3；电压幅值误差一般在0.001到0.01 pu。我建议初始设置取R = diag([1e-4, 1e-5])，角度方差的单位是rad^2，电压方差的单位是pu^2，这样噪声量级比较符合实际。

P0的初始化也容易踩坑。有些同学把P0设得特别小，觉得初始状态很准，结果前几步卡尔曼增益过小，估计值半天追不上真实轨迹。合理做法是给P0一个适中偏大的值，比如对角元素取1e-4到1e-3量级，让滤波器在前几步有足够的修正能力。

3.3 EKF主循环代码实现与雅可比处理

先写一个通用的数值雅可比函数，这样EKF主循环就不需要手推复杂的偏导公式了：

function J = numJacobian(f, x, varargin) % 中心差分数值雅可比矩阵 % f: 函数句柄，x: 状态向量，varargin: 其他参数 n = numel(x); fx = f(x, varargin{:}); m = numel(fx); J = zeros(m, n); h = 1e-6; % 扰动步长 for j = 1:n xp = x; xm = x; xp(j) = xp(j) + h; xm(j) = xm(j) - h; J(:, j) = (f(xp, varargin{:}) - f(xm, varargin{:})) / (2*h); end end

状态方程和量测方程写成独立函数，便于复用：

function dx = sysFun(x, u) % 三阶发电机模型状态方程 % x = [delta; dw; Eq] del = x(1); dw = x(2); Eq = x(3); Id = (Eq - Vinf*cos(del)) / xSum; Iq = Vinf*sin(del) / xSum; Pe = Eq*Vinf*sin(del) / xSum; dx = zeros(3,1); dx(1) = wb * dw; dx(2) = (Pm - Pe - D*dw) / M; dx(3) = (Efd - Eq - (xd - xdp)*Id) / Tdo; end function z = mesFun(x) % 量测方程：功角和机端电压幅值 del = x(1); Eq = x(3); Id = (Eq - Vinf*cos(del)) / xSum; Iq = Vinf*sin(del) / xSum; Vd = xq * Iq; Vq = Eq - xdp * Id; Vt = sqrt(Vd^2 + Vq^2); z = [del; Vt]; end

注意状态方程里我用的是全局变量或者外部参数，实际写代码时最好把参数定义在工作区或者用结构体传入，避免硬编码。

EKF主循环代码：

n = 3; m = 2; x_ekf = zeros(n, N); P = P0; x_ekf(:,1) = x0; for k = 1:N-1 % 状态预测：用Euler积分先跑通，后续可换RK4 x_pred = x_ekf(:,k) + dt * sysFun(x_ekf(:,k)); % 预测协方差 F = numJacobian(@sysFun, x_ekf(:,k)); P_pred = F * P * F' + Q; % 量测预测与雅可比 z_pred = mesFun(x_pred); H = numJacobian(@mesFun, x_pred); % 卡尔曼增益 S = H * P_pred * H' + R; K = P_pred * H' / S; % 更新 x_ekf(:,k+1) = x_pred + K * (z_meas(:,k+1) - z_pred); P = (eye(n) - K*H) * P_pred; end

这里有几个细节值得说明。第一，Euler积分在0.01s步长下勉强能用，但如果你发现结果有轻微振荡，不要急着怀疑滤波器，先把状态预测换成四阶Runge-Kutta，很多怪问题会自己消失。第二，雅可比矩阵用的是数值差分，效率不是最高，但对于三状态系统完全够用；上了多机系统再考虑解析式。第三，量测更新用的z_meas(:,k+1)是当前时刻的量测，不要错用成上一时刻。

3.4 UKF主循环代码实现与Sigma点生成

UKF的实现比EKF多一个Sigma点生成和权重计算的步骤：

% 参数设置 alpha = 0.01; 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) = lambda / (n + lambda) + (1 - alpha^2 + beta); for i = 2:2*n+1 Wm(i) = 1 / (2*(n + lambda)); Wc(i) = 1 / (2*(n + lambda)); end x_ukf = zeros(n, N); x_ukf(:,1) = x0; P = P0; for k = 1:N-1 % 1. 生成Sigma点 % 协方差矩阵正定修正 P_chol = P + 1e-12 * eye(n); S = chol(P_chol, 'lower'); chi = zeros(n, 2*n+1); chi(:,1) = x_ukf(:,k); for i = 1:n chi(:,i+1) = x_ukf(:,k) + sqrt(n + lambda) * S(:,i); chi(:,i+n+1) = x_ukf(:,k) - sqrt(n + lambda) * S(:,i); end % 2. Sigma点通过状态方程传播 chi_pred = zeros(n, 2*n+1); for i = 1:2*n+1 chi_pred(:,i) = chi(:,i) + dt * sysFun(chi(:,i)); end % 3. 加权合成预测均值和协方差 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 d = chi_pred(:,i) - x_pred; P_pred = P_pred + Wc(i) * (d * d'); end % 4. 量测更新：所有Sigma点通过量测方程 Z_pred = zeros(m, 2*n+1); for i = 1:2*n+1 Z_pred(:,i) = mesFun(chi_pred(:,i)); end z_mean = zeros(m,1); for i = 1:2*n+1 z_mean = z_mean + Wm(i) * Z_pred(:,i); end Pzz = R; % 量测协方差 Pxz = zeros(n,m); % 交叉协方差 for i = 1:2*n+1 dz = Z_pred(:,i) - z_mean; dx = chi_pred(:,i) - x_pred; Pzz = Pzz + Wc(i) * (dz * dz'); Pxz = Pxz + Wc(i) * (dx * dz'); end % 5. 增益与状态更新 K = Pxz / Pzz; x_ukf(:,k+1) = x_pred + K * (z_meas(:,k+1) - z_mean); P = P_pred - K * Pzz * K'; end

UKF比EKF多出来的计算量主要就在这两层循环上，但换来的是不需要任何求导操作，而且对非线性函数的逼近精度高得多。如果你发现chol(P_chol)仍然报错，说明协方差矩阵的状态很差，这时候把修正项从1e-12调到1e-6，基本能解决问题。

3.5 Q、R、P0的整定方法

滤波器的参数整定，本质上是给“模型预测”和“量测修正”分配信任度。Q大意味着模型不准、要更相信量测；R大意味着量测噪声大、要更相信模型。这个比例关系直接决定了滤波器的动态响应特性。

Q的初始值可以按状态量在正常运行时的动态范围来估计。比如功角的变化幅度可能在0.01到0.1 rad量级，那么对应的过程噪声方差可以取(0.01^2)到(0.1^2)之间的值，也就是1e-4到1e-2。角速度偏差和暂态电动势也类似。我常用的做法是：

Q = diag([1e-4, 1e-4, 1e-4]);

R的初始值按PMU量测精度取：

R = diag([1e-4, 1e-5]);

P0取一个相对适中的值，比如：

P0 = diag([1e-3, 1e-3, 1e-3]);

调参的总体原则是：先固定R，把Q从1e-6逐渐扫描到1e-2，观察估计误差的RMSE随Q的变化；再反过来固定Q扫描R。不要指望一次成功，多跑几组参数，你会看到明显的规律。后面第4节我会详细讲怎么从曲线形态反推参数问题。

4. 仿真结果分析与调参实战

4.1 扰动场景设置与结果观察

为了对比EKF和UKF在动态过程中的表现，我设置了一个典型的扰动算例：在t = 1s时刻，机械功率Pm从0.8 pu阶跃增加到0.84 pu，持续时间0.2s后恢复。这个扰动的物理含义相当于发电机输入功率的短时波动，幅度不大但足够让功角出现明显摆动。

仿真时长取5s，步长0.01s，总共500个量测点。真实轨迹用ode45生成，量测加噪声以后再分别输入EKF和UKF。跑完以后把功角、角速度偏差、暂态电动势三个状态的估计曲线和真实值画在同一张图上，比较直观。

从我的实测结果看，稳态阶段EKF和UKF都能跟上真实轨迹，误差差距不大；但扰动发生的瞬间，EKF的功角估计会出现一个相对明显的尖峰误差，且回到平稳状态的时间更长。UKF的估计曲线始终贴着真实值，尤其在功角变化率最大的那几个时刻，优势非常明显。这个现象和理论完全吻合：扰动越强，系统非线性越突出，EKF的一阶线性化误差就被放大得越厉害。

4.2 评价指标：RMSE与最大误差

单看曲线容易主观，还得用数值指标说话。我一般用均方根误差（RMSE）和最大绝对误差（MAE）两个指标：

rmse = sqrt(mean((x_est - x_true).^2, 2)); max_err = max(abs(x_est - x_true), [], 2);

RMSE反映整体跟踪精度，MAE反映最坏情况下的偏离程度。由于量测噪声是随机生成的，单次仿真的结果有随机性，我建议做20次蒙特卡洛仿真，每次重新生成噪声，然后统计RMSE的均值和标准差。这样得出来的结论才有说服力，审稿人也不会挑毛病。

我跑20次的平均结果显示：UKF的功角RMSE比EKF小约30%到50%，最大误差小得更多；在角速度偏差上，UKF的优势更大，因为角速度是由功角差分演变来的，状态方程非线性更敏感。如果你跑出来两者几乎一样，可以先检查是不是扰动设置太小，或者Q、R参数处于特殊比例上。

4.3 调参方向判断：从曲线形态反推问题

调参的时候，曲线会给你很多信息。我是这样判断的：

如果滤波曲线明显滞后于真实轨迹，在阶跃扰动后要好几拍才跟上去，通常是Q偏小或者R偏大，滤波器太信任模型、不敢采信量测。解决办法是把Q稍微调大，或者把R调小。

如果滤波曲线抖得很厉害，噪声完全没被平滑掉，甚至比量测噪声还剧烈，那就是R偏小或者Q偏大，滤波器把量测噪声当成了真实状态变化。解决办法是把R调大，让滤波器更信任模型。

如果估计值和真实值之间有一个稳定的偏置，扰动前后都存在，那不是滤波器的问题，多半是系统模型参数有偏差，或者量测方程写错了。这种系统性偏差靠调Q、R是救不回来的。

还有一个我很推荐的检查手段：打印“预测残差”，也就是量测值与量测预测值之差。如果滤波器调得合适，残差序列应该围绕0波动，且其协方差统计量接近R矩阵。如果残差明显有偏，或者方差远远大于R，说明模型和量测之间有一方出了问题。

我实际经验中还有一个性价比很高的做法：先做一个参数敏感性扫描，把Q的每个对角元素从1e-6到1e-2按对数间隔取5到6个值，跑完整的蒙特卡洛仿真，画出RMSE热力图。这个过程能让整个系统的行为规律一目了然，比手动盲调高效得多。

5. 常见问题与排错速查

5.1 高频故障与解决方案表

跑这套代码的时候，我遇到过的、以及帮别人排查过的问题集中在下面几类：

其中量测缺失这个场景值得多说一句。动态状态估计的优势之一就是能在量测短时缺失时维持输出，实现方法非常简单：把更新步的卡尔曼增益强制设为0，或者直接跳过更新部分，让状态和协方差用预测值。恢复量测后，滤波器会自动通过增益恢复正常更新，不需要额外逻辑。这个能力在PMU通信丢帧的情况下非常实用。

5.2 几条提高效率的实操建议

第一，代码结构上把状态方程、量测方程、数值雅可比、EKF主循环、UKF主循环都封装成独立函数，这样换系统模型的时候只需要改sysFun和mesFun，滤波器核心代码一行都不用动。我在从单机模型换到三机九节点系统的时候，只改了模型函数和参数，主循环完全复用，省了大量时间。

第二，调试阶段先用小步长、短时长、无扰动场景跑通，再逐步加扰动、加噪声。如果你一开始就直接上强扰动高噪声，出了问题都不知道该怀疑滤波器还是怀疑代码。

第三，运行蒙特卡洛仿真的时候，把噪声种子设为可配置的。你可以固定一个种子复现实验，也可以随机化种子做统计评估。Matlab里用rng(seed)控制即可。

第四，如果打算把动态状态估计作为论文研究点，建议在单机模型跑通以后，尽快过渡到IEEE标准测试系统，比如三机九节点或者IEEE 39节点系统。算法部分不需要改，但状态方程要从单机扩展到多机，量测矩阵维度会变大，需要考虑稀疏存储和更高效的雅可比计算。多机系统的状态量之间通过网络方程耦合，量测方程的维度也会显著增加，不过这已经超出本文范围，感兴趣的可以沿着这个方向继续深入。

还有一个细节，绘制结果图的时候，把量测噪声点也画出来，用浅灰色散点表示，可以直观看出滤波器的平滑效果。很多同学只画估计值和真值，忽略了量测噪声，结果图看起来过于“完美”，反而不真实。

最后再分享一个技巧：调完参数以后，把最终使用的Q、R、P0、alpha、beta、kappa这几个关键参数和对应的RMSE指标一起保存下来。这个习惯看似简单，实际帮了我大忙——同一个系统换扰动场景或者换量测配置时，我不需要从头盲调，直接在已有的参数网格上插值就行。做工程和做研究，能保存中间过程的人，永远比从头再来的人快。