去年冬天我拿到一篇关于扩展卡尔曼滤波在配电网故障测距中应用的论文本以为三个晚上就能把仿真图复现出来结果卡在状态方程推导上就耗了快两天。回头看卡住我的不是EKF本身的公式而是论文里大量“众所周知”的细节——量测方程里的故障电流怎么取、状态初值怎么给、仿真用多大的采样窗、滤波协方差怎么调这些论文往往只给结论不给路径。这篇文章不打算把论文复述一遍而是把从读完论文到跑通结果的全过程拆开讲模型怎么建、矩阵怎么写、数据怎么造、滤波为什么发散、最后怎么收敛。适合正在做毕业设计的电力方向研究生、需要评估算法落地效果的继电保护工程师以及所有想把仿真类论文真正跑起来的同行。1. 配电网故障测距为什么绕不开EKF1.1 测距问题的本质从故障量反推两个未知量配电网最常见的是单相接地故障故障后我们能在变电站测量端拿到三相电压和电流波形但真正想知道的是两个量故障点到测量端的距离以及故障点的过渡电阻。前者是测距结果后者虽然不是测距目标但它和距离混在一个非线性方程里不联合估计就没法把距离分离出来。这个问题的数学结构是故障距离d和过渡电阻Rf同时出现在故障回路电压方程里而且它们和电流、电压之间是乘积、相位耦合的非线性关系。传统阻抗法把故障电阻近似忽略或者单独估计在配电网这种多分支、负荷复杂、故障电阻不确定的场景下误差很容易被放大。EKF的价值在于它把“估计距离”和“估计电阻”放到同一个状态空间里用递推方式不断修正两个状态而不是先算一个阻抗再查表。1.2 EKF相比最小二乘和启发式算法的优势如果把故障后一小段窗口的数据拿来做批量最小二乘也能估计d和Rf但有个现实问题配电网故障暂态过程中电压电流相量并不是平稳的批量方法对数据窗口的起止时刻很敏感。EKF是逐拍递推的每个采样点或每个滑动相量点都能更新状态天然适合跟踪故障发生后的过渡变化。另外EKF的实现成本很低。和粒子滤波、无迹卡尔曼滤波相比它不需要采样粒子也不需要生成sigma点只需要一个状态方程、一个量测方程和一次雅可比矩阵求导。对论文复现来说这是巨大的优势——你把雅可比推对了后面的迭代逻辑基本就是标准卡尔曼公式代码量非常小调试复杂度也低。粒子滤波当然更鲁棒但配电网故障测距的论文里EKF仍然是出现频率最高的选项不是没有理由的。提示EKF适合的是“非线性程度中等、模型结构明确”的问题。如果量测方程极度非线性、多峰严重EKF会被线性化误差牵着走那就得考虑无迹卡尔曼或者粒子滤波了。复现之前先判断这一点能省很多调参时间。2. 复现基础先把状态方程和量测方程“焊死”2.1 状态量选取与基本假设我复现时把状态向量取成两个量(x [d, R_f]^T)其中d是故障距离单位kmRf是过渡电阻单位Ω。为什么只取这两个因为系统模型的状态转移很简单故障的物理参数在一个短时间窗内基本不变化所以状态方程可以写成(x_{k1} x_k w_k)也就是恒等转移加上一个很小的过程噪声。这里w_k的协方差矩阵Q用来吸收故障参数缓慢变化的部分。如果故障电弧电阻随时间漂移Q适当给大一点就能让滤波跟着走如果认为参数完全恒定Q就给得很小。这里必须注意一个前置假设故障类型已知且是单相接地。很多论文会先做故障选线、选相再进入EKF测距模块。复现时我建议先从最简单的A相单相接地开始跑通后再扩展到其他故障类型。2.2 一个可以直接落地的复数量测方程故障回路方程的复相量形式可以写成(\Delta \dot U_m d \cdot Z_l \cdot \Delta \dot I_m R_f \cdot \dot I_f)式中 (\Delta \dot U_m) 是测量端电压故障分量(\Delta \dot I_m) 是测量端电流故障分量(Z_l) 是线路单位长度正序阻抗(\dot I_f) 是故障支路电流。这个式子用文字描述就是测量端的故障电压由两部分构成一部分是沿线阻抗压降另一部分是过渡电阻上的压降。问题在于 (\dot I_f) 不容易直接测量。不同论文的处理方式不一样有的把它当成与 (\Delta \dot I_m) 成比例的电流有的把它的相角也放进状态向量。我复现入门用的简化假设是(\dot I_f \approx \Delta \dot I_m)也就是认为故障电流约等于测量端的电流故障分量。这个假设在单侧电源供电或者对侧电源较弱的配电网中成立得比较好复现时能很快跑通。但如果你的配电网模型是双端强电源这个假设会带来明显误差我后面会讲怎么扩展。把近似代入后量测方程变为(\Delta \dot U_m (d \cdot Z_l R_f) \cdot \Delta \dot I_m)2.3 从复数方程到实虚部展开代码里没法直接处理复数量测标准做法是把复数方程拆成实部和虚部两个实数方程。令(Z_l R_l jX_l)(\Delta \dot I_m I_x jI_y)(\dot I_f I_{fx} jI_{fy})(\Delta \dot U_m U_x jU_y)展开后得到两个量测方程(U_x d(R_l I_x - X_l I_y) R_f I_{fx})(U_y d(R_l I_y X_l I_x) R_f I_{fy})也就是说量测向量 (z [U_x, U_y]^T)状态向量 (x [d, R_f]^T)量测函数h(x)就是上面两个右式。这一步是整个复现的关键很多代码跑不出来不是因为EKF写错而是这里复数拆分的符号不对比如把 (-X_l I_y) 写成 (X_l I_y)滤波很快就会发散。提示不同论文对 (\dot I_f) 的定义会有差异一定要先去读目标论文里故障电流的说明。有些论文会把对侧电流也加进来有些会把故障电流相位单独作为一个状态量。复现的第一原则是论文假设什么你就实现什么不要轻易用自己以为“更合理”的模型替代。3. EKF递推与矩阵实现代码能跑起来的版本3.1 预测-更新循环的矩阵形式EKF的循环结构对于懂卡尔曼滤波的人来说不陌生但配电网测距这个场景有个特殊点状态转移矩阵F就是单位阵所以预测步简化为(x_{pred} x_k)(P_{pred} P_k Q)量测更新则是(K P_{pred} H^T (H P_{pred} H^T R)^{-1})(x_{k1} x_{pred} K(z - h(x_{pred})))(P_{k1} (I - K H) P_{pred})其中H是量测函数h对状态x的雅可比矩阵。在2.3节的模型下雅可比可以直接手推(H \begin{bmatrix} R_l I_x - X_l I_y I_{fx} \ R_l I_y X_l I_x I_{fy} \end{bmatrix})注意这个H矩阵每个时刻都在变因为Ix、Iy、Ifx、Ify都是当前时刻的相量值。这意味着我们不能像线性时不变系统那样离线算好增益K而是每个时刻都要重新计算。3.2 Python实现核心片段假设你已经通过滑动DFT得到了一个相量序列每个时刻的电压、电流实虚部分别存在U_seq、I_seq里故障电流用If_seq表示。核心循环如下import numpy as np # 线路参数有名值或标幺值均可示例为有名值 R_l 0.245 # 单位长度电阻 Ω/km X_l 0.355 # 单位长度电抗 Ω/km # 状态初值距离给线路中点过渡电阻给一个常见值 x np.array([5.0, 10.0]) # [d_km, Rf_ohm] P np.eye(2) * 1.0 # 过程噪声和量测噪声协方差 Q np.diag([1e-4, 1e-4]) R np.diag([1e-3, 1e-3]) # 运行EKF for k in range(len(U_seq)): Ux, Uy U_seq[k] # 当前电压相量实部、虚部故障分量 Ix, Iy I_seq[k] # 当前电流相量实部、虚部故障分量 Ifx, Ify If_seq[k] # 当前故障电流实部、虚部 # 预测步状态转移为恒等 x_pred x.copy() P_pred P Q d, Rf x_pred # 量测预测 h np.array([ d * (R_l * Ix - X_l * Iy) Rf * Ifx, d * (R_l * Iy X_l * Ix) Rf * Ify ]) # 雅可比矩阵 H np.array([ [R_l * Ix - X_l * Iy, Ifx], [R_l * Iy X_l * Ix, Ify] ]) # 更新步 z np.array([Ux, Uy]) y z - h S H P_pred H.T R K P_pred H.T np.linalg.inv(S) x x_pred K y P (np.eye(2) - K H) P_pred # 记录 d 和 Rf 的估计值 # d_est.append(x[0]); Rf_est.append(x[1])这段代码大约40行是完整可运行的EKF核心。我在复现时会把这段逻辑封装成一个类方便反复换数据测试。3.3 协方差初值、Q与R的调整口诀新手最容易卡死在Q和R的调整上。我的经验是先做量纲处理再做量级调节。如果直接拿有名值跑电压量级是10^4电流量级是10^2R矩阵对角元如果给成10^-3相当于告诉滤波“我对量测非常有信心”但实际传感器和DFT相位提取的误差远不止这个数滤波就会被模型误差带着走。建议把电压、电流先归一到标幺值再跑EKF。基值可以选额定相电压峰值和额定负荷电流峰值线路阻抗也除以对应阻抗基值。归一化之后状态量d和Rf都在0.01到几的量级Q和R的对角元从0.01、0.0001这样开始试马上就能找到合理范围。R的物理含义是量测误差主要由DFT相位提取误差决定给到电压基值的0.1%~1%比较合理。Q的物理含义是状态本身的漂移速度配电网故障参数短时间窗内很稳Q不要给太大否则协方差永远不会收敛估计结果会一直抖。我的初始值固定套路是P0取单位阵Q取对角线1e-4R取对角线1e-3然后看新息序列的统计特性微调。4. 仿真数据生成论文复现最花时间的环节4.1 用仿真模型造故障波形实测数据不好找论文复现几乎都是自己生成故障波形。我在Simulink里搭了一个10kV中性点经消弧线圈接地的配电网模型单端电源、一条10km线路、末端带2MW恒阻抗负荷故障点设在离测量端2km、5km、8km三个位置。核心参数如下表参数取值额定电压10kV线路长度10km单位长度电阻0.245 Ω/km单位长度电抗0.355 Ω/km故障类型A相单相接地故障电阻1Ω / 10Ω / 50Ω故障起始时间0.1s采样率10kHz采样率的选择有点讲究。基波相量测距用的是50Hz分量理论上一周期采20点就够但我实测发现采样率太低时DFT输出相位抖动很大EKF的收敛曲线会很毛糙。10kHz对应每周期200个采样点既能保证DFT精度又不至于数据量太大跑太慢。4.2 基波相量和故障分量的提取从原始波形到EKF能用的量测z中间隔着一步滑动DFT。我用的窗口是整一个工频周期即20ms每来一个新采样点窗口滑动一次计算当前窗口的基波相量(X(m) \frac{2}{N}\sum_{n0}^{N-1} x(n) e^{-j2\pi n/N})其中N200。这样得到的是随时间变化的相量序列每个时刻k对应一个U相量和一个I相量。故障分量的做法是取故障前0.04s到0.08s的平均基波相量作为参考值然后用故障后的每个相量减去这个参考值。为什么要扣故障前量因为负荷电流在线路阻抗上的压降也在测量电压里不扣掉的话量测方程里的“故障分量”含义就对不上了估计出的距离会整体偏差一段。4.3 数据对齐和起始点的门道还有一个容易翻车的细节故障发生后滑动DFT窗口内同时包含故障前和故障后的数据这时候DFT输出的是一个过渡值既不是故障前相量也不是故障后相量。如果直接从故障时刻开始跑EKF前10ms到20ms的量测数据本身就不可信滤波自然会乱跳。我的处理是从故障发生时刻往后推一个完整工频周期再加2ms也就是故障后22ms才开始启动EKF。这样保证窗内数据全部来自故障后的稳态过程。很多论文复现时省略了这个说明直接导致跑出来的距离曲线一开始有个巨大尖峰然后又慢慢爬回真值附近。这不是算法问题是数据起始点没选对。5. 复现踩坑调试从发散到收敛的完整链路5.1 现象一估计距离直接发散到乱跳我第一次跑EKF时d的估计值在前几个点就冲到上百km完全失控。当时的排查过程是这样的先打印每一拍的新息y的模值发现y从10^-3量级直接涨到10^2量级这说明量测预测h(x_pred)和实际量测z对不上问题大概率不在卡尔曼更新而在模型本身。关掉更新步把K矩阵强制设为0只让量测预测h跟着初值走发现h算出来的电压比实际电压大好几倍。这时候怀疑是不是线路阻抗的单位错了——Z_l必须是单位长度阻抗如果误用了全线路总阻抗h就会线性放大。检查后又发现复数拆分的符号问题电抗项前面漏了负号。修掉之后新息模值降到了合理范围滤波器开始收敛。这个经验是EKF发散的时候先别调Q和R先用“关更新看预测”的技巧定位是模型错了还是协方差给错了。5.2 现象二不发散但收敛到错误值第二类问题更隐蔽滤波器没发散d的估计值稳定在某个值附近但和真实故障距离对不上。我遇到过稳定收敛到比真实距离偏大约15%的情况。排查方向是先看故障分量扣干净没有。我当时用的是全量电压电流而不是故障分量负荷电流在故障前后变化会引入一个固定偏差这个偏差直接变成d的系统误差。换成故障分量序列后收敛值立刻回到了真实距离附近。另外一个常见来源是初值给得太离谱。EKF本质是局部线性化如果状态初值和真值偏离太远雅可比矩阵在初值点的线性化可能把迭代方向引向局部极值。我的经验是d的初值给线路中点Rf初值给5到10Ω不要给0也不要给几百欧否则一旦掉进局部极值后续测量更新很难拉回来。5.3 现象三波形曲线和论文对不上代码跑通、估计距离也基本对了但画出来的收敛曲线和论文里的图轮廓对不上。这种问题多半出在故障起始角上。故障初相角决定了故障电压突变时刻在工频周期里的位置直接影响暂态过程和基波相量的初始相位。论文里如果用的是电压峰值时刻初相角90度起故障而你用的是电压过零点两条曲线自然对不上。我建议复现时把故障起始角统一设成90度附近这样暂态分量最强基波相量提取反而更稳定。另外故障电阻使用不同的值也会改变收敛速度论文可能在50Ω工况下画图你用10Ω去跑收敛时间差一倍都不奇怪。先跑几组不同的Rf和故障距离确认规律一致再对齐具体数值这样才算真正理解了论文的结论而不是只得到一张长得像的图。5.4 把新息序列当调试窗口调试EKF最该盯着看的不是状态估计曲线而是新息序列 (y_k z_k - h(x_{pred}))。理论上如果模型正确、协方差合理新息应该是零均值、白噪声特性的序列。我每次跑完一组数据都会画一下新息曲线和它的直方图新息表现判断均值明显偏离0量测方程有系统偏差检查故障分量和阻抗参数方差远大于R模型误差大或Q过小方差远小于RQ过大状态被过度信任可能出现“假收敛”有明显的趋势变化状态真值在漂移Q给大一些或检查数据起始点这个方法比盯着d的曲线可靠得多。状态曲线看着收敛有可能是协方差假象新息才是滤波器和数据之间的真实契约。6. 从复现别人的论文到改进自己的算法6.1 把简化假设放开加入故障电流相角2.2节的简化假设 (\dot I_f \approx \Delta \dot I_m) 在双端强电源系统下不成立。改进方向是把故障电流的幅值和相角作为额外的状态量或者只把相角θ放进去状态向量变成(x [d, R_f, \theta]^T)量测方程变成(\Delta \dot U_m d \cdot Z_l \cdot \Delta \dot I_m R_f \cdot |\dot I_f| e^{j\theta})此时量测方程对θ的偏导不再为零雅可比矩阵变成3x2结构滤波精度会更高但代价是状态维数增加初值更难给。我在做双端仿真时用过这个扩展故障电阻估计抗噪能力明显提升。这个扩展是论文复现之后很自然的演进方向也是把“复现”变成“研究”的起点。6.2 论文复现的通用套路模型假设是灵魂复现过配电网测距论文后再看其他仿真类复现任务比如一些用COMSOL做激光熔覆仿真的复现工作流程其实是相通的第一步把论文里的数学假设列出来第二步把参数表整理成配置文件第三步先搭最小可运行原型第四步逐步逼近论文的仿真条件。很多人复现失败是因为一开始就追求完全一致结果被细节淹没。先跑通最小原型再一层层加条件才是最快的路径。6.3 后续还能做的方向EKF测距本身已经很成熟但配电网场景下可改进的方向仍然不少。比如EKF的Q、R矩阵固定值不一定最优可以用Sage-Husa自适应滤波在线估计噪声协方差比如配电网含分布式电源后双向潮流让单端量测方程不再可靠可以尝试双端量测融合再比如高阻接地时故障分量太小EKF容易抖动可以结合小波变换先做特征提取再进滤波器。这些方向不需要更换整个算法框架都是在现有EKF结构上的增量修改。先把基础复现跑通剩下的问题就都是明确的问题而不是盲目的试错。我个人在实际复现中最大的体会是复现一篇论文更像替作者做一次代码审查。你以为卡住你的会是卡尔曼滤波那套公式实际上真正花时间的往往是论文里被一句话带过的模型假设、参数取值和数据预处理。把这些细节啃下来之后你会突然发现原来论文里每个量测方程的系数都不是随手写的背后都对应着物理场景里的一个实际约束。这就是复现能带来但单纯读论文得不到的东西。