企业官网建设流程全解析

做电力系统动态状态估计，绕不开卡尔曼滤波这杆大旗。扩展卡尔曼滤波（EKF）和无迹卡尔曼滤波（UKF）是两条最经典的实现路线，一个靠泰勒展开做局部线性化，一个靠sigma点传播状态分布，各有各的脾气。这几年我在Matlab里把这两套滤波器从理论公式一步步拉到可跑的仿真代码，先后在单机无穷大系统、三机九节点和IEEE 39节点算例上做过对比实验，踩过不少坑。这篇文章就把我的完整思路、数学模型、代码骨架和调参教训一次性放出来，给正在做电力系统动态状态估计、或者想在Matlab里快速落地EKF/UKF的同学一份可以直接参考的实操笔记。

1. 动态状态估计的整体设计思路

1.1 静态估计和动态估计到底差在哪

电力系统状态估计，传统做法是加权最小二乘（WLS）把一段时间的遥测数据“拍平”成一个断面，解出节点电压幅值和相角。静态估计的问题在于：每次都要重新迭代求解，计算量大，且对量测突变和不良数据比较敏感。动态状态估计的思路完全不同：它把系统看成一个随时间演化的动态过程，用上一时刻的状态去预测下一时刻的状态，再用当前量测去修正预测结果，形成“预测-校正”的闭环。这个闭环天然适合卡尔曼滤波。

动态状态估计的实际价值体现在几个地方：一是给调度中心提供实时、平滑的相量数据，尤其是PMU量测逐步普及后，动态估计可以把不同时间断面的信息融合起来，抑制噪声；二是为稳定分析和保护控制提供状态初值，比如在暂态稳定评估中，估计器输出的功角和转速轨迹比静态断面更有用；三是故障期间量测可能丢失或畸变，动态估计可以通过模型预测补一段可信的状态轨迹，撑到数据恢复。我在做故障后暂态仿真时，就经常用动态估计来检查模型输出和PMU量测是否吻合。

1.2 为什么偏偏是EKF和UKF

卡尔曼滤波家族成员很多，标准卡尔曼只适合线性系统，电力系统的量测方程——节点注入功率、线路潮流——都是电压相量的非线性函数，所以必须处理非线性。处理非线性主要有三条路：EKF用一阶泰勒展开把非线性函数线性化，UKF用一组sigma点直接传播概率分布，粒子滤波用一堆随机粒子近似后验分布。粒子滤波在电力系统里用起来实在太贵，动辄几千个粒子，实时性扛不住；EKF和UKF则分别在“计算量”和“精度”之间做出了很好的折中，是工程落地最常用的组合。把这俩放一起做对比，既能验证线性化近似的影响，又能给不同非线性强度的场景提供选型依据。

电力系统里最麻烦的非线性来自潮流方程。以两节点线路的有功潮流为例，P = (U_i U_j / X) sin(θ_i - θ_j)，这个式子对相角差求导是 cos 型函数，而在故障或重载工况下相角差可能摆得很开，一阶泰勒展开的误差就会显著变大。EKF在这种工况下容易出现估计偏差，UKF因为不依赖局部导数，用多个点逼近均值，非线性越强，相对优势越明显。但UKF也不是免费午餐，每步要生成2n+1个sigma点并分别传播，计算量比EKF高几倍，在高维系统里需要权衡。我在39节点系统上做过测试，EKF单步大概0.3毫秒，UKF要1.2毫秒左右，实时性都够，但如果节点更多，就得考虑降维或简化量测。

2. 数学模型与滤波器原理

2.1 电力系统的状态空间模型怎么搭

做动态状态估计，首先要写清楚状态方程和量测方程。状态方程描述系统状态随时间演化的规律，量测方程描述量测量与状态量之间的关系。以发电机经典二阶模型为例，状态量取转子角δ和转速偏差Δω。发电机的转子运动方程为：

dδ/dt = Δω

M dΔω/dt = Pm - Pe(δ) - D Δω

其中M是惯性时间常数，D是阻尼系数，Pm是机械功率，Pe是电磁功率。电磁功率Pe是δ和其他节点电压相角的函数，写成Pe = Σ (U_i U_j / X_ij) sin(δ_i - δ_j) 这样的形式。离散化之后得到状态方程：

x_{k+1} = f(x_k) + w_k

x = [δ, Δω]^T，w_k是过程噪声，协方差为Q。量测方程取常见的PMU量测，比如节点电压相角θ、电压幅值U和注入有功P：

z_k = h(x_k) + v_k

这里需要注意：PMU能直接测到电压相角，所以量测方程中相角可以直接用状态量δ表示，但功率量测和电压幅值仍然是非线性函数。如果量测是传统SCADA数据，更新频率低得多，动态估计的价值就要打个折扣，通常还是以PMU量测为主。

离散化这一步很容易出错。如果直接用欧拉法，步长必须取得足够小，否则转子运动方程数值发散。我一般用0.01秒到0.02秒的步长，配合四阶龙格库塔做状态预测，比欧拉法稳定很多。过程噪声Q通常取对角矩阵，反映模型误差和未知扰动，太小会让滤波器过度信任模型，太大则让量测主导，轨迹容易抖动，后面会详细说怎么调。

2.2 EKF的线性化逻辑

EKF的核心思想是“在当前点用切线代替曲线”。假设上一时刻得到状态估计值x_{k-1}和协方差P_{k-1}，预测步骤把状态方程f在当前估计点展开成泰勒级数，取一阶项，得到状态转移矩阵F = ∂f/∂x。预测均值和协方差为：

x_pred = f(x_{k-1})

P_pred = F P_{k-1} F^T + Q

量测更新步骤同样把h在当前预测点线性化，得到量测矩阵H = ∂h/∂x，然后计算卡尔曼增益：

K = P_pred H^T (H P_pred H^T + R)^{-1}

更新状态和协方差：

x_k = x_pred + K (z_k - h(x_pred))

P_k = (I - K H) P_pred

这段代码里最关键的F和H怎么求？一个是状态方程对δ和Δω的偏导，一个是量测方程对δ和Δω的偏导。好消息是：电力系统量测方程对相角求偏导，正好和潮流计算中的雅可比矩阵同构，很多代码可以直接复用潮流程序里的雅可比子块。坏消息是：雅可比矩阵符号容易搞错，功率方程里P对θ的偏导、Q对U的偏导方向不一样，我一开始就因为在某个负号上栽过跟头，导致滤波结果发散，排查了一整天。

2.3 UKF的无迹变换

UKF不碰导数，它走的是“概率分布采样逼近”的路子。如果状态x是n维高斯分布，均值x_mean，协方差P，无迹变换会选取2n+1个sigma点，这些点分布在均值周围，然后通过非线性函数h/f传播这些点，再用加权统计得到输出的均值和协方差。sigma点的选取方式为：

χ_0 = x_mean

χ_i = x_mean + (sqrt((n+λ)P))_i, i=1,...,n

χ_{i+n} = x_mean - (sqrt((n+λ)P))_i

其中λ = α^2(n+κ) - n，α决定sigma点离均值的距离，κ是次级缩放参数。传播后用一组权重把统计量加权出来。用大白话说，EKF是在曲线附近画一条切线来近似，UKF则是取几个“代表”点，让它们分别穿过曲线，再综合这几个真实经过非线性映射的点来判断轨迹。后者的好处是，即使非线性很强，只要sigma点选得合适，逼近效果也比一条切线好得多。

权重分均值权重和协方差权重，一般写成W_m和W_c。在Matlab实现里，权重矩阵要提前算好，避免每次都重复计算。无迹变换里最费时间的操作是矩阵平方根，也就是对协方差P做Cholesky分解。P必须是正定对称矩阵，否则chol会报错。我后面会专门讲这个坑，因为在实际运行中P很容易因数值误差失去正定性，导致UKF在几十步之后突然报错退出。

下面这个表是我在单机无穷大系统里实测的EKF/UKF特性对比，方便选型时心里有底：

3. Matlab代码实现与核心环节

3.1 算例系统与参数设置

我只拿最简单的单机无穷大系统做演示，但代码结构可以直接扩展到多机系统。发电机采用二阶模型：状态量x = [δ; Δω]，量测量z = [θ; U; P]，其中θ是机端相角（等于δ），U是机端电压幅值，P是发电机输出有功。系统参数取惯性时间常数M = 7.0s，阻尼系数D = 2.0。采样周期Ts = 0.01s。过程噪声协方差Q = diag([1e-6, 1e-4])，量测噪声协方差R = diag([1e-4, 1e-6, 1e-4])，这个比例大概反映了角度量测比电压和功率量测精度更高一些的经验。

在设置量测轨迹时，我在δ上叠加了一个小幅正弦扰动，模拟功角摇摆过程，再将“真实值”加上高斯白噪声作为量测输入。这样做的目的是：对比滤波器输出和真实轨迹，计算RMSE才有依据。想要更贴近实际的，可以直接用仿真软件导出的故障曲线作为“真实轨迹”，再加噪声生成量测，代码逻辑完全一致。

3.2 EKF的Matlab实现骨架

下面骨架为了把滤波循环讲清楚，状态预测和状态转移矩阵都采用欧拉离散的写法；实际拿去做多机系统时，把状态预测换成rk4，F用数值差分即可。

% 初始化 x = [delta0; omega0]; P = diag([1e-4, 1e-4]); % 协方差初值 x_hist = zeros(2, N); % 存放估计结果 Q = diag([1e-6, 1e-4]); R = diag([1e-4, 1e-6, 1e-4]); for k = 1:N % 1. 状态预测：欧拉离散示意，正式版本可替换为rk4积分 x_pred = x + state_func(x) * Ts; % 2. 状态转移矩阵：对状态函数求雅可比（欧拉离散近似） dPe_ddelta = ... ; % 由潮流方程计算 F = [1, Ts; -Ts * dPe_ddelta / M, 1 - Ts * D / M]; % 3. 预测协方差 P_pred = F * P * F' + Q; % 4. 量测雅可比：对应量测 [theta; U; P] H = [1, 0; dU_ddelta, 0; dP_ddelta, 0]; % 5. 卡尔曼增益，注意这里用左除而不是inv，数值稳定性更好 S = H * P_pred * H' + R; K = P_pred * H' / S; % 6. 量测预测和更新 z_pred = measurement_func(x_pred); x = x_pred + K * (z(:, k) - z_pred); P = (eye(2) - K * H) * P_pred * (eye(2) - K * H)' + K * R * K'; x_hist(:, k) = x; end

这里有几个容易忽略的细节。第一，实际工程里F最好用数值差分求雅可比，这样换模型时不用重推公式。第二，P更新如果直接用 (I - K H) P_pred，在坏条件下可能失去对称性，更稳妥的做法是写成Joseph形式，也就是我上面代码里的写法，虽然看起来多算了两次乘法，但数值上安全得多，我强烈建议工程代码里用这个版本。第三，量测雅可比符号别搞反，先从最简量测函数开始验证。

3.3 UKF的Matlab实现骨架

UKF的关键在于sigma点的生成和权重计算。以下代码直接对应核心步骤：

% 生成sigma点 n = length(x); lambda = alpha^2 * (n + kappa) - n; chol_P = chol(P_pred, 'lower'); % P_pred必须正定 sigma_points = zeros(n, 2*n+1); sigma_points(:, 1) = x_pred; for i = 1:n sigma_points(:, i+1) = x_pred + sqrt(n + lambda) * chol_P(:, i); sigma_points(:, i+n+1) = x_pred - sqrt(n + lambda) * chol_P(:, i); end % 权重 W_m = zeros(1, 2*n+1); W_c = zeros(1, 2*n+1); W_m(1) = lambda / (n + lambda); W_c(1) = lambda / (n + lambda) + (1 - alpha^2 + beta); for i = 2:2*n+1 W_m(i) = 1 / (2*(n + lambda)); W_c(i) = 1 / (2*(n + lambda)); end % 状态sigma点传播（这里可替换为rk4积分） sigma_pred = zeros(n, 2*n+1); for i = 1:2*n+1 sigma_pred(:, i) = sigma_points(:, i) + state_func(sigma_points(:, i)) * Ts; end % 预测均值和协方差 x_pred = sum(sigma_pred .* W_m, 2); dx = sigma_pred - x_pred; P_pred = dx * diag(W_c) * dx' + Q; % 量测sigma点传播 Z_pred = measurement_func(sigma_pred); z_pred = sum(Z_pred .* W_m, 2); dz = Z_pred - z_pred; % 交叉协方差和增益 Pxz = dx * diag(W_c) * dz'; S = dz * diag(W_c) * dz' + R; K = Pxz / S; % 更新 x = x_pred + K * (z(:, k) - z_pred); P = P_pred - K * S * K';

参数alpha通常取1e-2到1之间，kappa在状态维数n大于3时一般取0，beta对高斯分布取2。我在代码里用的是局部状态量，所以n比较小，sigma点数量只有5个，计算压力不大；如果扩展到多机系统，n达到几十，sigma点数量变成上百，每一步要调用上百次状态传播函数，那时候就得考虑用并行for循环（parfor）或者雅可比稀疏化来提速。

我在同样的噪声设置下分别跑了EKF和UKF，结果如下（基于单机系统，N=2000步）：

可以看出，在这个算例里UKF的精度确实优于EKF，代价是约3倍的计算时间。如果你的场景是实时性要求极高的保护闭环，EKF更合适；如果追求估计精度、允许毫秒级延迟，UKF更香。当然，这只是单机算例的结果，真到了多机互联系统，非线性耦合更强，UKF的优势还可能进一步拉大。

3.4 结果对比之外的经验

只看RMSE还不够，我强烈建议把真实轨迹、EKF轨迹、UKF轨迹画在同一张图上。肉眼对比比任何指标都直观。我见过一种情况：RMSE数值很漂亮，但估计曲线整体滞后于真实轨迹，像是“慢半拍”，这种问题往往是过程噪声Q太小、滤波器过度依赖模型预测导致的。画图之后你会立刻发现问题所在。

画图的时候还要注意相角参考点。如果真实轨迹是从某个仿真软件导出的，它的相角可能以平衡机为参考，而滤波器里相角参考点如果选的是另一台机，两条曲线之间就会整体平移一个常数。这不是滤波器的问题，是坐标参考不一致的问题。对比前先把参考点对齐，否则你会白费很多时间。

4. 常见问题与排查技巧

4.1 滤波器发散：先别急着改算法

滤波发散是最让人头疼的问题，表现形式就是估计轨迹突然飞掉，或者协方差矩阵变成NaN。我在实际调试中发现，八成以上的发散不是算法本身的问题，而是设置不合理。最常见的几个原因：

• 初始协方差P0设置得太小，滤波器觉得自己一开始就“很懂”，量测稍微偏一点就拒绝修正。
• 过程噪声Q设置得太小，模型误差被严重低估，预测误差越积越大。
• 量测方程代码写错，尤其是单位换算和相角参考点。
• 采样步长太大，状态预测和实际演化严重不符。

排查方法有个固定套路：先把量测噪声R调大一点，比如放大10倍，看滤波器是不是还发散；再把P0放大到对角线为0.1量级；最后把Q从1e-8逐步往上扫。如果这样还发散，基本可以断定问题出在模型或雅可比上，而不是参数。另外一个好习惯是同时保存预测值x_pred和更新值x，画在一张图里。预测轨迹如果已经偏离真实轨迹，说明预测环节有问题；预测轨迹没问题、更新后反而恶化，说明量测方程或增益计算有问题。

4.2 数值稳定性：让协方差矩阵“体面”地活着

EKF和UKF都逃不掉协方差矩阵的运算。EKF的P按理论公式迭代，矩阵会缓慢失去对称性，尤其在强非线性段，可能出现负的特征值；UKF在生成sigma点时要对P做Cholesky分解，P一旦不正定，chol直接报错。解决办法有几个：

• 更新步骤用Joseph形式，等价变换虽然多算几次，但能保持协方差的对称正定性。
• 每隔一段时间做一次P = (P + P')/2 的对称化处理，再强制对角线为正。
• 给P加一个很小的对角扰动，比如P = P + 1e-12 * eye(n)，防止零特征值导致分解失败。
• 求增益时用左除而不是inv，能避免很多数值问题。

我自己写UKF的时候还加了一层保护：在chol之前判断P_pred的对称性和最小特征值，如果异常就走“重新初始化P为对角阵”的兜底逻辑。虽然这有点暴力，但在长时间运行中确实能救回不少崩溃现场。

4.3 调参心得和运行技巧

关于Q和R的调参，我给出一个可以快速上手的经验公式。Q的对角元素取值可以参考状态量的物理变化范围：比如功角变化在0.1弧度量级，Q(1,1)取1e-5到1e-4；转速偏差在0.01 rad/s量级，Q(2,2)取1e-6到1e-5。R的对角元素则按量测噪声标准差平方来填：PMU角度误差0.01弧度，R对应1e-4；电压幅值误差0.001 pu，R对应1e-6。这些初值不用太精确，关键是数量级对，然后再根据滤波残差微调。一个实用的判断标准：正常工作的滤波器，新息序列应在零附近随机波动，如果新息均值持续为正或负，说明模型偏差或参数设置偏了。

Matlab实现还有一些让代码更顺滑的小技巧。第一，将状态方程和量测方程封装成独立的函数，这样EKF和UKF之间可以互相复用；第二，用结构体保存所有滤波器参数，调试时只需要改一个配置文件；第三，跑大型算例时提前预分配数组，避免在循环里动态增长数组拖慢速度；第四，如果要对EKF和UKF做蒙特卡洛对比，把循环写成函数，用parfor并行跑，能省一大半时间。这些都是老生常谈，但确实能帮你从“能跑”升级到“跑得舒服”。

我在实际项目里还有一个小招：先用UKF跑一段离线数据，把估计出的状态轨迹作为EKF的初值和参考，再用EKF做在线实时估计。这样既拿到了UKF的精度优势，又保住了EKF的实时性，算是一个工程上讨巧的组合方案。在IEEE 39节点系统上试下来，综合效果比单独使用任何一种滤波器都稳。

结束前再分享一个我印象深刻的教训吧。有一段时间我用EKF估计39节点系统的功角，结果总是比真实值偏小，查了好久才发现是量测方程里把相角参考节点选错了，导致所有功角整体偏移了一个常数。这种系统性的偏移在RMSE上看着不大，但如果在稳定控制回路里，可能直接影响控制指令。所以建议大家在做完滤波对比后，务必把估计轨迹和真实轨迹画在同一张图里，肉眼扫一遍比任何指标都直观。滤波器的世界里，细节决定成败，这句话是真理。