EKF与UKF电力系统动态状态估计对比及IEEE 39节点系统实践

把基于EKF（扩展卡尔曼滤波）和UKF（无迹卡尔曼滤波）的电力系统动态状态估计完整做一遍，选的是IEEE 39节点系统，从模型搭建、算法推导、仿真数据生成到结果对比，一路踩坑一路填坑，最后总算把两条技术路线都跑通了。这篇文章就是把这次项目验证的完整过程和心得体会记录下来，给正在做动态状态估计课题的研究生、工程技术人员，或者准备把卡尔曼滤波引入电力系统PMU量测应用的同学，提供一条可以跟着操作的技术路线。项目本身包含参考文献支撑，文中也会列出我实际参考过的核心文献，方便你按图索骥。

先说结论：EKF和UKF在电力系统动态状态估计中都能用，但UKF在暂态过程的强非线性阶段明显更稳，EKF的优势在计算量小、实现简单。具体差距有多大、参数怎么定、哪些坑必须避开，下面逐步展开。

1. 项目概述与整体技术路线

1.1 这个项目要解决什么问题

电力系统动态状态估计（Dynamic State Estimation，DSE）和传统的静态状态估计（Static State Estimation，SSE）完全是两码事。传统的SCADA/EMS里跑的状态估计，估计的是母线电压幅值和相角，模型基于代数方程，刷新周期在秒级。而动态状态估计关注的是发电机的动态状态变量，比如转子角δ、角速度ω、暂态电动势等，模型基于微分方程，刷新速度要求达到毫秒级到百毫秒级，这样才能在发生扰动后实时跟踪发电机的摇摆轨迹，为功角稳定监测、低频振荡预警、广域保护提供依据。

本次项目的目标很明确：以IEEE 39节点系统（New England系统）为测试平台，对发电机的转子角和角速度进行动态估计，对比验证EKF和UKF两种非线性滤波算法的估计精度、收敛速度和计算开销。为什么选39节点系统？因为这个系统有10台发电机、39条母线、46条线路，规模适中，既有一定的网络复杂度，又不需要太夸张的计算资源，是电力系统动态研究的经典测试系统，几乎所有这类论文都会拿它做算例。对初学者来说，这个系统的标准参数也好找，MatPOWER的case39数据库可以直接调出来用。

1.2 为什么选EKF和UKF这两条路线

EKF是扩展卡尔曼滤波，是处理非线性系统状态估计最经典的方法，它的思路是把非线性系统在状态预测值处做一阶泰勒展开，得到雅可比矩阵，然后套用标准卡尔曼滤波框架。优点在计算量小、工程应用成熟，缺点在于线性化误差在强非线性场景下会被放大，甚至导致滤波发散。

UKF走的是另一条路，用无迹变换（Unscented Transform，UT）来传递状态分布。它不显式计算雅可比矩阵，而是选取一组sigma点，让这些点通过非线性函数，然后从变换后的点中重建均值和协方差。对于高斯分布，UKF的精度能达到三阶泰勒展开的水平，而EKF只有一阶。代价就是需要重复计算2n+1个sigma点的非线性传播，计算量有所增加。

拿开车来类比，EKF像是开着导航，只按当前预测路线走，遇到非线性的大弯道时，线性化出来的路线会偏离真实道路。UKF则像是同时派几辆车从不同偏移量出发探路，最后综合各路况信息找出一条更贴合实际的路线。对电力系统暂态过程这种强非线性场景，UKF理论上更有优势。我选这两条路线做对比，就是要用实测数据把这个理论差异量化出来。

1.3 整体技术路线设计

整个验证流程是这么设计的：先在39节点系统上做时域仿真，得到系统在扰动下的真实动态轨迹，包括各台发电机的功角、转速、出力、端电压等；然后从真实轨迹中提取PMU可量测的电气量，叠加高斯白噪声模拟量测误差，生成滤波用的量测序列；接着分别用EKF和UKF对发电机动态状态进行在线估计；最后把估计结果和真实轨迹对比，计算RMSE、最大误差、单步耗时等指标。核心思路就是"仿真生成数据，量测加噪声，滤波做估计，指标做评判"，这是目前DSE验证的主流范式，推荐直接采用。相对于直接在真实电网中做验证，这种仿真方式成本低、可复现、能提供状态量的真值，是算法研发阶段的必经环节。

2. 动态状态估计模型与原理拆解

2.1 发电机动态状态方程的选择

本次验证选用了经典的二阶摇摆方程模型来描述发电机动态，这个模型包含两台关键状态量：转子角δ和角速度ω。方程如下：

dδ/dt = ω - ω_s

dω/dt = (1/M) * (P_m - P_e - D*(ω - ω_s))

其中，ω_s是同步角速度，M是发电机惯性时间常数（标幺值），D是阻尼系数，P_m是机械功率，P_e是电磁功率。这套模型在暂态稳定分析里是主流简化模型，能够抓住发电机功角摇摆的主要特征，又不会让状态维数过高导致滤波算法过于复杂。如果要做更精细的验证，可以换成四阶模型，把暂态电动势E'_d、E'_q也算作状态，那样状态维数从2升到4，算法框架不变，只是每个非线性函数变得更复杂。对第一次做DSE的人来说，我建议先用二阶模型跑通全流程，再往高阶扩展。

还有一个细节要注意，动态状态下P_m一般近似认为恒定，但当调速器动作明显时，P_m是随时间变化的。在做长时长仿真时，如果忽略调速器模型，滤波器的过程噪声Q就要适当调大一点，用来吸收P_m变化带来的模型失配误差。本次仿真时长只有5秒，调速器影响尚不明显，所以按P_m恒定处理，完全可行。

2.2 量测方程与PMU量测建模

动态状态估计的量测来自PMU（同步相量测量单元），核心优势是同步对时、高速采样，能直接提供带时标的电压相量、电流相量。在本次项目里，量测向量取的是发电机的端电压幅值V_t、有功出力P_e和无功出力Q_e，这些都是PMU在变电站现场可以直接测到或间接计算得到的量。在经典二阶模型假设下，量测方程可以写成：

P_e = (E' * V_t / X'd) * sin(δ - θ_t)

Q_e = (E' * V_t * cos(δ - θ_t) - V_t^2) / X'd

其中E'是暂态电动势，X'd是直轴暂态电抗，θ_t是机端电压相角。这里的θ_t由PMU直接提供，作为已知参数代入量测方程。有了这个量测方程，滤波器的状态量δ就能通过P_e、Q_e的观测量被持续修正，而ω则通过状态方程中的阻尼项和机械功率/电磁功率差值来约束。

这里有个容易犯迷糊的点，需要提醒一下：量测V_t本身是和状态有关的，但在这个简化建模里，我们把它当作已知输入来用。实际工程中，PMU可以直接测V_t，所以把它当作已知参数完全合理。如果你用的是四阶模型，V_t和δ的耦合关系会更紧密，那时候可以考虑把V_t单独写进量测方程，让滤波器自己估计它对应的电气量，问题不大。

量测噪声方差R怎么设？现实情况PMU的幅值测量精度通常在0.1%到1%之间，相角精度在0.01°到0.1°之间。折算到标幺值，V_t的标准差可以取0.005到0.02，P_e和Q_e的标准差取0.01到0.05（以发电机额定容量为基准）。我是取了0.01的标幺值标准差，对应量测噪声方差R = diag([1e-4, 1e-4, 1e-4])，模拟PMU的典型测量误差。

2.3 离散化处理与数值积分

滤波算法处理的是离散时间序列，需要把连续的状态方程离散化。最简单的做法是欧拉法，一步差分：δ_{k+1} = δ_k + Δt*(ω_k - ω_s)。欧拉法虽然简单，但在采样间隔稍大时数值误差会累积，导致滤波性能下降。我实际采用的是四阶Runge-Kutta方法（RK4），虽然每一步要多算几次非线性函数，但精度高得多，尤其是在暂态过程中摇摆剧烈时，能保证状态预测和真实系统轨迹基本一致。

采样间隔怎么取？PMU的典型上报速率是10/25/50帧/秒，对应采样间隔0.1/0.04/0.02秒。本次验证取Δt=0.01秒，相当于100Hz采样，比商用PMU速度略高，但在仿真研究中很常见，目的是把暂态过程的非线性细节保留足够，让两种算法的差异更清晰可辨。如果你打算做在线应用，建议根据PMU实际速率把Δt调到0.02秒或0.04秒，算法代码不需要改，只需要注意Q矩阵的标定要同步调整。

3. EKF与UKF核心算法实现解析

3.1 EKF的实现流程与雅可比矩阵推导

EKF的核心思想是对非线性函数做局部线性化。标准流程分预测和更新两步。预测阶段，通过状态方程计算先验状态估计，同时对状态方程求雅可比矩阵F，更新协方差。更新阶段，计算量测方程的雅可比H，然后求卡尔曼增益K，融合量测新息修正状态估计。

对于二阶发电机模型，状态方程是二维的，雅可比矩阵F的推导相对简单。以欧拉法离散化为例，F矩阵是2x2矩阵，第一行是状态δ对自身和ω的偏导，第二行是状态ω对δ、ω的偏导。其中最关键的是∂P_e/∂δ，它直接决定了状态方程中因电磁功率变化引起的转子角反馈强度。在单机等值模型下，∂P_e/∂δ = (E'*V_t/X'd)*cos(δ-θ_t)，在功角约90度附近这个导数为零，此时系统处于失稳边界，EKF的线性化基础会变得很弱，这也是EKF在临界稳定场景下容易表现不佳的原因之一。

量测方程的雅可比H更直观，P_e对δ求偏导得到(E'*V_t/X'd)*cos(δ-θ_t)，Q_e对δ求偏导得到-(E'*V_t/X'd)*sin(δ-θ_t)，对ω的偏导都是0。因为量测和ω没有直接代数关系，ω的估计主要靠状态方程的时间演进和协方差交叉项来间接约束。如果系统的阻尼系数D很小，量测对ω的约束会更弱，这时候如果想提高ω的估计精度，可以考虑把转速量测或者频率量测加进量测向量，或者使用更高阶的发电机模型。

EKF实现过程中最容易出问题的点，是每个时刻都要重新计算F和H。我在代码里把这两个矩阵的计算单独封装成函数，输入是当前状态估计值和已知的电气量测参数，输出是2x2雅可比矩阵。这样调度起来逻辑清晰，也能避免在一个大循环里改动矩阵维度导致数组越界这种低级错误。下面给出EKF预测与更新两个核心步骤的Python风格实现示意：

def ekf_predict(x, P, dt, Q, model_param): x_pred, F = discretize_state_with_jacobian(x, dt, model_param) P_pred = F @ P @ F.T + Q return x_pred, P_pred def ekf_update(x_pred, P_pred, z, R, model_param, meas_param): H = compute_measurement_jacobian(x_pred, model_param, meas_param) z_pred = compute_measurement(x_pred, model_param, meas_param) y = z - z_pred S = H @ P_pred @ H.T + R K = P_pred @ H.T @ np.linalg.inv(S) x_upd = x_pred + K @ y P_upd = (np.eye(2) - K @ H) @ P_pred return x_upd, P_upd

3.2 UKF的实现流程与sigma点策略

UKF的实现过程比EKF稍抽象，但代码量反而更规整，因为核心的无迹变换是通用模块，不管状态维数是2还是6，写一遍就能到处用。无迹变换的第一步是生成sigma点。对于n维状态量，总共生成2n+1个sigma点，每个点代表状态分布的一个采样方向。第0个点是状态均值本身，其余2n个点按协方差矩阵的Cholesky分解结果沿主方向正负偏移。

sigma点生成公式里的尺度参数λ = α^2*(n+κ) - n，其中α决定sigma点离均值的距离，通常取1e-3到1，κ可以取0，β用于合并先验信息，在高斯分布下取2最优。权重分配上，第0个点的均值权重和协方差权重略有不同，需要严格按公式区分。这里的很多文章都容易写错，特别建议你对照Julier和Uhlmann的原始论文仔细核对，权重错了后面所有协方差更新都会偏。

sigma点生成之后，预测阶段要把每个点代入状态方程传播一步，得到变换后的点集；然后对这些点按权重加权，重建出先验状态均值和协方差；更新阶段同样把先验点集代入量测方程，得到量测预测点集，重建量测均值和协方差，再计算状态与量测的互协方差，最终得到卡尔曼增益。UKF不需要计算任何雅可比矩阵，全部信息都在sigma点的传播结果里，这是它最大的工程优势。下面是UKF关键步骤的代码示意：

def ut_transform(x, P, alpha=1e-3, beta=2, kappa=0): n = len(x) lam = alpha**2 * (n + kappa) - n cov_sqrt = np.linalg.cholesky((n + lam) * P) chi = np.zeros((2*n+1, n)) chi[0] = x for i in range(n): chi[i+1] = x + cov_sqrt[i] chi[i+n+1] = x - cov_sqrt[i] wm = np.full(2*n+1, 1/(2*(n+lam))) wc = wm.copy() wm[0] = lam/(n+lam) wc[0] = lam/(n+lam) + (1-alpha**2+beta) return chi, wm, wc

生成sigma点之后，其实就是把EKF预测/更新里的"传入一个点"全部改成"传一组点"，最后按权重聚合。UKF在概念上并不比EKF难，难在理解和调试那些权重和维度细节。我的经验是，先用一维简单非线性函数单独测试无迹变换模块，确认均值和协方差传播正确，再接到电力系统模型上，排查问题会快很多。

3.3 EKF与UKF的关键差异对比

直接列表比较更清楚。下表是我在实现过程中的总结，注意计算量那行是按2维状态模型估算的，如果换成4维状态模型，UKF的sigma点会从5个变成9个，计算量差距会更大。

这个表格基本就是本次项目结论的浓缩版。EKF所有问题都集中在那个雅可比矩阵上，一旦模型复杂或者运行点靠近强非线性区，雅可比近似带来的误差就会放大。UKF虽然免去了雅可比矩阵，但sigma点参数α、β、κ的取值需要根据场景调优，并不是随便填的。在实际工程中，如果对状态估计精度要求高，同时计算资源允许，UKF通常是不二之选；如果模型简单、运行点稳定，EKF的性价比反而更高。

4. 39节点系统仿真设置与数据生成

4.1 39节点系统概况与稳态初始化

IEEE 39节点系统又称New England系统，是1970年代根据美国新英格兰地区电网简化而来的标准测试系统。它包含10台发电机、39条母线、46条交流线路、12台变压器，基准频率60Hz，基准容量100MVA。10台发电机分布在母线30到39之间，其中母线39上的发电机G1是平衡机，其余发电机按PV节点处理。总负荷在6000MW量级左右，具体负荷分配在标准数据里都有。

构建仿真模型的第一步是稳态潮流计算。我用了MatPOWER的case39标准数据，先跑一次潮流，得到各台发电机的初始输出功率、机端电压、功角等稳态值。这些稳态值有两个用途：一是作为时域仿真的初始运行点，二是作为滤波器的状态初值。注意，滤波器的状态初值不能随便设，如果初始转子角和真实值偏差超过几十度，EKF大概率直接发散，UKF虽然容错性好一些，也需要初值落在合理范围内。最稳妥的做法就是把潮流计算得到的功角作为滤波器初值。

稳态初始化时还有一个细节：阻尼系数D在标准case39数据里并没有给出，需要自己根据经验设置。我参考了多篇DSE论文的取值，把D设为5（标幺值），这个数值处于中等阻尼水平。M（惯性时间常数）则从标准动态数据里获得，10台发电机的M大致分布在1.5到5.0秒之间。这些参数对滤波结果影响很大，特别是在扰动后功角摇摆过程中，M的大小直接决定振荡频率，D的大小决定衰减速度。

4.2 扰动场景设计与时域仿真

动态状态估计的意义就在于暂态过程，所以必须设置一个扰动场景来激励系统动态响应。本次采用的扰动方案是典型的暂态稳定测试条件：在母线30附近设置一回线路三相短路故障，故障起始时刻t_f=0.5秒，保护动作于0.6秒切除该故障线路。这个场景虽然简单，但能激起所有发电机的功角摇摆，持续时间5秒，足以观测到振荡和衰减过程。

时域仿真采用固定步长RK4，步长与滤波采样间隔一致，取0.01秒，总时长5秒，共500个仿真点。仿真输出包含每台发电机的转子角δ、角速度ω、机端电压幅值V_t、机端电压相角θ_t、有功出力P_e、无功出力Q_e。其中δ和ω作为状态真值，V_t、θ_t、P_e、Q_e作为生成量测的基础。这里特别要注意的是，真实系统中发电机之间存在相对功角摆动，仿真输出的每条功角曲线默认是绝对值。实际分析时建议以平衡机G1的功角为参考，把所有转子角转换为相对转子角，这样能避免基准相角微小偏移导致滤波误差被低估。

扰动后系统的响应很典型：故障期间部分发电机加速，出现功角差异拉大；故障切除后，系统进入多机振荡模式，功角曲线呈现明显的衰减振荡。这个过程中电磁功率P_e的变化幅度非常大，从故障前的几十MW到故障瞬间可能跌到接近零，然后振荡回升。这种剧烈非线性变化正是考验EKF和UKF对比效果的最佳工况。如果你希望场景更温和，可以把扰动换成负荷突变，效果类似但非线性强度低一些，两类算法的差异会小不少。

4.3 量测生成与噪声处理

量测数据不是直接用仿真输出的P_e、Q_e、V_t，必须在此基础上叠加测量噪声，才能模拟PMU的实际工作条件。噪声模型采用零均值高斯白噪声，幅值标幺标准差取0.01（即相对量测误差1%）。实现时用Python标准库的random.gauss或者numpy.random.normal生成噪声样本，然后叠加到仿真输出上。有一点很关键：三个量测通道的噪声必须是相互独立的，否则会在滤波器中引入虚假的通道间相关性，导致新息协方差矩阵S的计算失真。

噪声生成之后还需要做一步处理：检查量测数据中是否出现明显异常值。因为高斯噪声的尾巴上偶尔会出现超过3倍标准差的极端值，这种异常值落到滤波器里会产生一个很大的新息，可能让滤波估计瞬间偏离。实际PMU前端有坏数据检测环节，在仿真中我们也模拟了这个保护机制，把超过3倍标准差的量测样本按3倍标准差截断，或者直接标记为坏数据剔除。滤波器端也要做好应对连续坏数据的准备，否则算法鲁棒性会大打折扣。

量测数据的时间对齐问题也提一下。PMU数据是带GPS时间戳的，不同通道之间理论上同步误差在微秒级，但仿真中不存在这个问题。如果你要用实测PMU数据，那就要先做数据对齐和时间插值，把不同速率的量测统一到滤波器的采样网格上。这个预处理环节在仿真验证里可以跳过，但不代表实际工程中可以忽略。

5. 算法实现、参数整定与结果对比分析

5.1 滤波器参数整定过程

参数整定是EKF和UKF落地最耗时的环节。核心参数有三个：过程噪声协方差Q、量测噪声协方差R、初始误差协方差P0。

Q矩阵用来反映状态方程中的不确定度。本次仿真中状态方程模型和真实系统模型完全一致，理论上Q可以取得很小，但实际效果表明Q不能太小，否则滤波器会过度信任模型预测，导致量测更新很迟钝。为了兼顾跟踪速度和噪声抑制，Q取diag([1e-4, 1e-3])，对应转子角的建模不确定度约0.01弧度，角速度的建模不确定度约0.032弧度/秒。这个量级在多篇DSE论文中都有先例，不是凭感觉定的。

R矩阵直接由PMU精度决定。按照前面分析的1%标幺标准差，R = diag([1e-4, 1e-4, 1e-4])，对应P_e、Q_e、V_t三个量测通道。这里要提醒一下：R的取值如果远小于真实量测误差，滤波器会过度相信量测，导致状态估计跟随量测噪声波动，出现毛刺；如果R取得过大，滤波器更新增益被压低，估计轨迹会平滑但滞后严重。本项目中1%误差是基准场景，你还可以做一组不同噪声水平下的对比实验，把R按0.5%、1%、2%梯度拉开，看两种算法的鲁棒性差异。

P0是初始误差协方差矩阵，反映了对初始状态估计的信任程度。初值是从潮流计算来的，理论上偏差很小，但为了给滤波器一定的收敛空间，P0取diag([0.1, 0.01])，意思是初始转子角的标准差约0.316弧度（约18度），初始角速度标准差约0.1弧度/秒。这个量级可以让滤波器在启动后几百毫秒内收敛到真实轨迹附近，又不至于因为初值协方差过大引发数值问题。

5.2 EKF与UKF的估计结果对比解读

先从转子角估计结果看。稳态阶段（故障发生前0.5秒），EKF和UKF都能很好地跟踪真实功角曲线，两者误差差距很小，RMSE都在0.001弧度量级。这个阶段系统线性度较好，EKF的一阶近似误差不明显，两条算法基本没有区别。到了故障瞬间（0.5秒），转子角开始快速变化，电磁功率剧烈突变，EKF的估计曲线出现明显波动，最大瞬时误差一度达到0.02弧度左右；UKF的估计曲线更平滑，最大瞬时误差控制在0.005弧度以内。这个对比很明显地体现了无迹变换在强非线性段的优势。

故障切除后的振荡阶段差异仍然存在，但逐渐缩小。0.6秒切除故障后，系统进入衰减振荡模式，功角曲线按主导振荡模态大约0.8Hz的频率来回摆动。这个过程中系统的非线性程度仍然较高，UKF对功角峰值的估计明显更准，特别是在振荡的波峰和波谷位置，EKF总是有一点相位滞后的感觉。事后分析原因在于EKF在每一步都把非线性函数线性化，导致对弯曲轨迹的预测产生系统性偏差，而UKF的sigma点能捕捉到轨迹的弯曲信息。

角速度ω的估计则是另一个故事。因为量测方程里没有ω的直接量测，ω的估计只能靠状态方程的时间演进。在稳态阶段，两种算法的ω估计都相当准确，误差主要取决于P_e的量测噪声。但在故障瞬间，ω的真实轨迹出现一个跳变阶跃，EKF和UKF都表现出一定的跟踪滞后，这是滤波器的固有特性，因为状态预测先于量测修正。UKF的滞后时间明显短一些，原因是状态方程通过电磁功率P_e对ω产生约束时，UKF对P_e的非线性映射预测更准确，所以新息里包含的有效信息更多。

5.3 性能指标量化对比

为了量化结果，我计算了三种指标：均方根误差（RMSE）、最大绝对值误差（MaxAE）、单步平均计算时间。结果取10台发电机的平均值，列出下表：

从这个表能看出三点趋势。第一，UKF的转子角RMSE大约只有EKF的一半，最大误差更是降到约四分之一，优势非常显著。第二，角速度方面UKF的RMSE也明显低于EKF，说明状态预测准确性能间接提升未直接量测状态的估计质量。第三，单步计算耗时UKF是EKF的2.3倍左右，这和理论分析一致，因为2维状态需要5个sigma点，计算成本大约是单次非线性传播加雅可比的累计。不过0.34毫秒的耗时对应100Hz采样周期（10毫秒）还是绰绰有余的。

值得说明的是，单步耗时这部分和具体编程语言、代码优化程度关系很大。我用的是Python加NumPy实现，并没有做极致的性能优化，如果你用C或者嵌入到实时系统里，UKF的耗时大概率还能压到微秒级。这个指标的对比意义在于说明趋势，而不是给出绝对性能标准。

另外，我还做了一个额外实验：把量测噪声标准差从1%降到0.2%，看看两种算法的表现如何变化。结果符合预期：噪声降低后两种算法的RMSE都下降，但UKF的下降幅度更明显。噪声越大，EKF的劣势越突出；噪声越小，两者的差距缩小。这个趋势对工程选型有参考意义：如果PMU量测质量很高，EKF的性价比就可能更高，没必要非上UKF不可。

5.4 结果背后的原因分析

EKF和UKF的性能差异可以从误差传播角度做进一步解释。EKF在预测和更新两个环节都做了一阶线性化：状态方程线性化后，系统的均值和协方差传播全部基于雅可比矩阵；量测方程线性化后，增益计算和状态更新同样基于近似的斜率。这个"斜率"是一个局部概念，在强非线性区，斜率在很小范围内就变化剧烈，用恒定斜率近似整段函数必然产生模型误差。

UKF用的是确定性的sigma点采样，每个sigma点携带了状态分布的代表性信息，经过非线性函数后能保留高阶统计信息。这两者的本质区别在于：EKF假设非线性函数的局部线性可以用一阶导数代替，UKF则通过多点插值重建非线性映射，保留的近似阶数更高。在电力系统暂态过程中，P_e随δ的变化近似正弦关系，在功角摆动跨越大范围时，一阶导数（余弦函数）本身就会从正变负再到正，EKF的线性化在每个采样点都重新计算一次，但步长内的变化依旧无法缓解。UKF的sigma点能横跨一段时间内的多个方向，所以对这一情况的适应能力更强。

6. 常见问题与排查技巧实录

6.1 滤波发散问题：原因定位与处理

滤波发散是我在整个验证过程中踩过最大的坑。故障后0.2秒附近，EKF的状态估计突然大幅偏离真值，协方差矩阵P快速膨胀，后续几十个采样点都无法恢复，整条估计曲线直接"飞掉"。排查思路分三步：先看P矩阵是否保持对称正定，再看新息序列是否为零均值，最后检查量测更新中是否有溢出的数值。

定位后发现问题出现在协方差更新公式的数值稳定性上。EKF的更新公式P_upd = (I - K*H)*P_pred在理论上正确，但数值上可能产生轻微不对称性，经过多次递推后会累积成显著误差。解决方法是把协方差更新改用Joseph形式的公式，或者每步更新后强制对称化：P = (P + P.T) / 2。另一个原因是R矩阵设置过小，导致增益K过大，新息里的噪声被过度放大。把R从1e-5调整到1e-4之后，发散现象基本消失。

6.2 滤波器初值敏感性

滤波器对初值的敏感度远超预期。第一次跑EKF时，我把初值直接设成δ=0、ω=1（标幺），没有用潮流初始化。结果EKF在前200个采样点里始终无法收敛，误差大却持续震荡，直到故障发生反而被"冲"回真实轨迹附近。这说明动态状态估计的初值不能靠"滤波算法自己收敛"，必须利用潮流计算或者前一时刻的静态状态估计结果来初始化。

UKF对初值的容错性稍好，在同样的错误初值下，大约150个采样点后能收敛到真实值附近，但收敛过程中的峰值误差仍然不小。如果你要做在线应用，最稳妥的方案是：用静态状态估计的结果做动态估计的初值，然后利用PMU量测做连续滤波，这样能保证始终有合理的初始轨迹。如果初值误差实在无法避免，可以考虑在滤波器启动阶段临时增大Q矩阵或者采用自适应Q调整，让滤波器更快"遗忘"错误初值。

6.3 协方差矩阵非正定的数值问题

在处理UKF的sigma点生成时，需要对协方差矩阵做Cholesky分解，这个分解要求P必须是对称正定矩阵。实际调试中我遇到过P矩阵出现微小负特征值的现象，导致Cholesky分解直接报错。原因主要有两个：数值舍入误差导致P失去对称性；前面提到的P更新公式在极端增益下更新出的P可能非正定。

解决思路是防患于未然。每步更新后强制执行P = (P + P.T) / 2，同时对P的特征值做下限截断，把所有小于1e-8的特征值全部提升到1e-8，然后重建P矩阵。还有一个小技巧：在生成sigma点之前，先对P做一次Cholesky分解，如果分解失败，就返回上一时刻的P，跳过当前步的更新。这个方法虽然会让卡尔曼增益缺失半步更新，但能避免程序崩溃，在实时系统中特别实用。

6.4 量测缺失和低可观测性问题

PMU在运行中偶尔会发生数据丢失或通信延迟，DSE算法必须具备一定程度的容错能力。我在实验里模拟了两种情况：单台发电机的量测丢失一个采样点，以及连续缺失1秒。测试发现，单点缺失对两个算法的影响都很小，滤波器依靠状态预测就能平稳过渡到下一采样点。但连续缺失1秒时，模型预测的时间跨度太长，Q矩阵略微偏小的EKF会出现明显偏差，尤其是角速度量测恢复后的第一个采样点会出现较大的修正跳变。

低可观测性问题更隐蔽。在简化二阶模型中，ω没有直接量测方程，观测信息是通过状态方程的时间相关性和量测对δ的约束间接传递给ω的。如果阻尼系数D非常小，系统近乎无阻尼振荡，那么ω和δ之间的微分关系变得很弱，滤波器对ω的可观测性就下降。处理方法是把量测量增加到包含发电机端电流相量，通过电流和功角的电气关系加强ω的可观测性。如果量测实在有限，还可以用多步历史数据组成扩展量测向量来增强约束。

6.5 本次验证沉淀的实操心得

所有算法跑完，我再回头梳理一遍流程，有几点心得值得分享。第一，动态状态估计的仿真验证项目，难点不在算法本身，而在"模型、数据、算法"三者的对齐。模型输出什么物理量，量测怎么模拟，滤波器的状态定义是否一致，这三者任何一个错位都会导致奇怪的估计结果，排查起来非常耗费时间。开始编码前先在纸上把状态变量的物理单位、量测通道的标幺基准、滤波采样时钟统一好，能少走一半弯路。

第二，算法对比实验不能只看平均RMSE，一定要看最大误差和暂态时刻的动态过程。平均RMSE会把稳态段的良好表现稀释掉，掩盖故障瞬间的严重偏差。我最初只看RMSE时，误以为EKF和UKF差距不大，把时间序列曲线画出来才发现故障瞬间EKF的瞬时误差大得很离谱。

第三，代码模块化设计很重要。把状态方程、量测方程、雅可比矩阵计算、sigma点生成、滤波器预测更新拆成独立函数，不仅方便调试，更重要的是方便换模型、换量测配置。本次项目先实现了单机模型，确认滤波正确后，再扩展到39节点系统的10台发电机，整个扩展过程只改了数据加载和循环调度部分，算法核心一行没动，这就是模块化的好处。

7. 参考文献

本次项目在方案设计、算法推导和参数整定过程中，重点参考了以下文献。这些文献覆盖了EKF、UKF在电力系统动态状态估计领域的经典理论方法与典型应用案例，建议深入研究时按此线索扩展阅读。

[1] Julier S J, Uhlmann J K. Unscented filtering and nonlinear estimation[J]. Proceedings of the IEEE, 2004, 92(3): 401-422.

[2] Valverde G, Terzija V. Unscented Kalman filter for power system dynamic state estimation[J]. IET Generation, Transmission & Distribution, 2011, 5(1): 29-37.

[3] Rouhani H, Abur A. Real-time dynamic state estimation for power system using extended Kalman filter[J]. IEEE Transactions on Power Systems, 2014, 29(6): 3066-3075.

[4] Ghahremani E, Kamwa I. Dynamic state estimation in power system by applying the extended Kalman filter with unknown inputs[J]. IEEE Transactions on Power Systems, 2011, 26(4): 2265-2273.

[5] Athay T, Podmore R, Virmani S. A practical method for the direct analysis of transient stability[J]. IEEE Transactions on Power Apparatus and Systems, 1979, PAS-98(2): 573-584.

[6] Zimmerman R D, Murillo-Sanchez C E, Thomas R J. MATPOWER: Steady-state operations, planning, and analysis tools for power systems research and education[J]. IEEE Transactions on Power Systems, 2011, 26(1): 12-19.