☰
从建模到调参:EKF配电网故障测距仿真复现全攻略
2026/10/5 7:41:11 网站建设 项目流程

去年冬天,我拿到一篇关于扩展卡尔曼滤波在配电网故障测距中应用的论文,本以为三个晚上就能把仿真图复现出来,结果卡在状态方程推导上就耗了快两天。回头看,卡住我的不是EKF本身的公式,而是论文里大量“众所周知”的细节——量测方程里的故障电流怎么取、状态初值怎么给、仿真用多大的采样窗、滤波协方差怎么调,这些论文往往只给结论,不给路径。这篇文章不打算把论文复述一遍,而是把从读完论文到跑通结果的全过程拆开讲:模型怎么建、矩阵怎么写、数据怎么造、滤波为什么发散、最后怎么收敛。适合正在做毕业设计的电力方向研究生、需要评估算法落地效果的继电保护工程师,以及所有想把仿真类论文真正跑起来的同行。

1. 配电网故障测距为什么绕不开EKF

1.1 测距问题的本质:从故障量反推两个未知量

配电网最常见的是单相接地故障,故障后我们能在变电站测量端拿到三相电压和电流波形,但真正想知道的是两个量:故障点到测量端的距离,以及故障点的过渡电阻。前者是测距结果,后者虽然不是测距目标,但它和距离混在一个非线性方程里,不联合估计就没法把距离分离出来。

这个问题的数学结构是:故障距离d和过渡电阻Rf同时出现在故障回路电压方程里,而且它们和电流、电压之间是乘积、相位耦合的非线性关系。传统阻抗法把故障电阻近似忽略或者单独估计,在配电网这种多分支、负荷复杂、故障电阻不确定的场景下误差很容易被放大。EKF的价值在于,它把“估计距离”和“估计电阻”放到同一个状态空间里,用递推方式不断修正两个状态,而不是先算一个阻抗再查表。

1.2 EKF相比最小二乘和启发式算法的优势

如果把故障后一小段窗口的数据拿来做批量最小二乘,也能估计d和Rf,但有个现实问题:配电网故障暂态过程中,电压电流相量并不是平稳的,批量方法对数据窗口的起止时刻很敏感。EKF是逐拍递推的,每个采样点或每个滑动相量点都能更新状态,天然适合跟踪故障发生后的过渡变化。

另外,EKF的实现成本很低。和粒子滤波、无迹卡尔曼滤波相比,它不需要采样粒子也不需要生成sigma点,只需要一个状态方程、一个量测方程和一次雅可比矩阵求导。对论文复现来说这是巨大的优势——你把雅可比推对了,后面的迭代逻辑基本就是标准卡尔曼公式,代码量非常小,调试复杂度也低。粒子滤波当然更鲁棒,但配电网故障测距的论文里EKF仍然是出现频率最高的选项,不是没有理由的。

提示:EKF适合的是“非线性程度中等、模型结构明确”的问题。如果量测方程极度非线性、多峰严重,EKF会被线性化误差牵着走,那就得考虑无迹卡尔曼或者粒子滤波了。复现之前先判断这一点,能省很多调参时间。

2. 复现基础:先把状态方程和量测方程“焊死”

2.1 状态量选取与基本假设

我复现时把状态向量取成两个量:

(x = [d, R_f]^T)

其中d是故障距离,单位km;Rf是过渡电阻,单位Ω。为什么只取这两个?因为系统模型的状态转移很简单:故障的物理参数在一个短时间窗内基本不变化,所以状态方程可以写成:

(x_{k+1} = 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_{k+1} = x_{pred} + K(z - h(x_{pred})))

(P_{k+1} = (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^2,R矩阵对角元如果给成10^-3,相当于告诉滤波“我对量测非常有信心”,但实际传感器和DFT相位提取的误差远不止这个数,滤波就会被模型误差带着走。

建议把电压、电流先归一到标幺值再跑EKF。基值可以选额定相电压峰值和额定负荷电流峰值,线路阻抗也除以对应阻抗基值。归一化之后,状态量d和Rf都在0.01到几的量级,Q和R的对角元从0.01、0.0001这样开始试,马上就能找到合理范围。

R的物理含义是量测误差,主要由DFT相位提取误差决定,给到电压基值的0.1%~1%比较合理。Q的物理含义是状态本身的漂移速度,配电网故障参数短时间窗内很稳,Q不要给太大,否则协方差永远不会收敛,估计结果会一直抖。我的初始值固定套路是:P0取单位阵,Q取对角线1e-4,R取对角线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_{n=0}^{N-1} x(n) e^{-j2\pi n/N})

其中N=200。这样得到的是随时间变化的相量序列,每个时刻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结构上的增量修改。先把基础复现跑通,剩下的问题就都是明确的问题,而不是盲目的试错。

我个人在实际复现中最大的体会是:复现一篇论文,更像替作者做一次代码审查。你以为卡住你的会是卡尔曼滤波那套公式,实际上真正花时间的往往是论文里被一句话带过的模型假设、参数取值和数据预处理。把这些细节啃下来之后,你会突然发现,原来论文里每个量测方程的系数都不是随手写的,背后都对应着物理场景里的一个实际约束。这就是复现能带来但单纯读论文得不到的东西。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询