简介这份资源围绕小电流接地系统单相接地故障选线这一电力系统保护难题展开面向从事继电保护与控制研究的研究人员、工程师及研究生。针对传统选线方法在参数变化、高阻故障与不同合闸角下准确率偏低的问题资源系统复现了论文提出的三种新原理基于工频与特征频率零序导纳相角差判据、综合低频与主频分量的改进Prony分解方法以及消弧线圈并联中值电阻系统电阻投切前后工频增量系数法并给出配套代码与逐步解释。包内为1个pdf文件约910KB内含论文复现分析、方法原理推导、Matlab仿真流程与等效电路建模说明可对照代码理解特征频率搜索、DCMP与APM指标构建及故障线路判定逻辑。已有71人学习适合希望掌握选线判据实现与仿真验证、为装置选型和技术方案制定提供参考的读者。1. 小电流接地选线为什么在参数变化下失灵从三种新判据说起配电网里中性点不接地或经消弧线圈接地的系统单相接地时故障电流常常只有几安到几十安现场叫小电流接地系统。传统选线靠零序电流比幅比相最头疼的是过渡电阻一高、故障合闸角一小、消弧线圈补偿度一调幅值和相位跟着漂选线装置就开始乱报。这个复现包围绕论文里的三条新判据展开工频与特征频率零序导纳相角差、改进 Prony 分解后的 DCMP 与 APM 双判据以及消弧线圈并联中值电阻投切前后工频增量系数 K。代码以 Python 仿真为主配套 Matlab 风格的信号生成适合做保护算法、选线装置整定和论文复现的人拆着看。2. 零序导纳相角差判据特征频率搜索与相角差整定2.1 判据的物理依据与适用范围小电流接地系统的零序网络里健全线路的零序导纳主要来自对地电容工频下阻抗角接近 -90°故障线路的零序电流要流过接地点过渡电阻和消弧线圈高频段还叠加线路电感与对地电容构成的谐振分量。论文抓住的点是同一条线路在工频和特征频率下分别算零序导纳相角故障线路两者的角差接近 180°健全线路接近 0°。这条判据不依赖零序电流的绝对幅值对高阻故障比传统比幅法稳。它的边界也很明确特征频率必须真的存在。特征频率一般在 150~250 Hz 范围具体落点跟线路长度、电缆比例、消弧线圈补偿度有关如果系统里电缆很短、过渡电阻又特别大这段高频能量可能被 CT 和滤波器吃掉。另外各条线路 CT 的变比、极性必须一致否则相角差本身没有可比性。观测量健全线路典型值故障线路典型值整定注意工频零序导纳相角接近 -90°容性受消弧线圈影响上移取稳定窗平均特征频率零序导纳相角与工频接近与工频相差接近 180°特征频率要能量占优相角差绝对值40°120°门槛可先取 120°2.2 S 变换特征频率搜索的代码实现与坑原复现里的 S 变换写法window np.exp(-2*np.pi**2*f**2*t**2)直接拿t np.arange(N)/fs构造高斯窗等价于把窗中心放在 t0结果是对信号起始段加权最重对稳态段几乎没有分辨力。S 变换的窗中心应该随每个采样点移动工程上要快速复现常见做法是用带通滤波加滑动高斯窗做近似先找能量最大的候选频点。import numpy as np from scipy.signal import firwin, lfilter, hilbert def search_characteristic_freq(signals, fs, f_min150, f_max250, step5): 在候选频段内搜索总能量最大的频率signals 形状 (n_lines, n_samples) best_f, best_e None, -np.inf for f in np.arange(f_min, f_max step, step): # 带通取候选频点附近 ±5 Hz避免全带噪声参与能量统计 b firwin(129, [f - 5, f 5], pass_zeroFalse, fsfs) total 0.0 for sig in signals: y lfilter(b, 1, sig) env np.abs(hilbert(y)) # 希尔伯特包络近似该频点瞬时能量 total np.sum(env ** 2) if total best_e: best_f, best_e f, total return best_f, best_e逻辑说明每个候选频率做一次带通用希尔伯特包络平方累积能量哪个频点总能量最大就选为特征频率。参数说明fs是采样率现场录波常见 4 kHz、10 kHz、12.8 kHzf_min/f_max按系统电容电流水平先粗估10 kV 配网一般在 150~250 Hz35 kV 或长电缆系统可能要往上挪到 300 Hz滤波器阶数 129 对应 10 kHz 采样约 12.9 ms 窗长低于 2 个工频周期够用但会带来约 6.5 ms 群延迟后续算相角时必须做延迟补偿。提示能量搜索不要用原始信号的 FFT 峰值直接代替FFT 峰值容易落在工频泄漏的旁瓣上。2.3 工频与特征频率相角差的提取与判据整定知道特征频率之后接下来对每条线路分别提取工频和特征频率相量。稳妥做法是整周期 DFT 或窄带滤波加希尔伯特。def band_phase(sig, fs, f0, bw5): 提取 f0 附近分量的解缠相位 b firwin(129, [f0 - bw, f0 bw], pass_zeroFalse, fsfs) y lfilter(b, 1, sig) z hilbert(y) return np.unwrap(np.angle(z)) def phase_diff_select(signals, fs, f_char, f050): n_lines signals.shape[0] diff np.zeros(n_lines) for i, sig in enumerate(signals): th_f band_phase(sig, fs, f0) th_c band_phase(sig, fs, f_char) k len(sig) // 3 # 跳过滤波器启动段 d np.angle(np.exp(1j * (th_f[k:] - th_c[k:]))).mean() diff[i] abs(d) return diff, int(np.argmax(diff))逻辑说明band_phase先带通再希尔伯特np.angle(exp(1j*Δ))把相位差折到 (-π, π]避免直接相减出现 2π 跳变。参数说明bw5是半带宽频率成分密集时可缩到 2~3 Hz但滤波器阶数要跟着加k用来跳过启动段。判据整定建议最大角差超过 120° 才输出选线结果最大与次大差值小于 20° 时报“选线不确定”不要硬指一条。参数建议值调整方向特征频率搜索步长5 Hz频带拥挤时降到 2 Hz带通半带宽2~5 Hz噪声大取大分辨率要求高取小相角差门槛120°现场误报多就抬到 140°最大/次大间隔20°小于该值降级为不确定3. 改进 Prony 分解与 DCMP/APM 双判据实现3.1 为什么用 Prony 而不是直接 FFTFFT 受窗长限制10 kHz 采样、0.1 s 窗才 10 Hz 分辨率低频段和衰减直流分不开Prony 把信号建模成若干复指数之和能给出频率、阻尼、幅值和相位短窗更合适。小电流接地故障的零序电流里同时有工频残流、低频衰减分量和高频分量低合闸角时低频分量比高频更突出论文用 DCMP 抓这个特征高合闸角时主频分量的幅相差异更明显用 APM。这两种场景如果只用 FFT 的幅值比就容易被频率泄漏拉平。Prony 的代价是对噪声和模型阶数敏感阶数给大了会冒出虚假极点频率和阻尼都是假的。工程上不会直接套原始 Prony常见做法是 SVD 降阶加总体最小二乘或者先做带通再分解。3.2 Hankel 矩阵与极点求解的可运行写法原复现代码里V1 V[:, :n]、V2 V[:, 1:n1]、A pinv(V1[:-1]) V2[:-1]这一段下标逻辑是拧的V是右奇异向量构成的矩阵按那个写法求出的 A 既不是方阵也没有明确的极点含义。更接近标准 SVD-Prony 的写法如下。import numpy as np def prony_decompose(x, fs, max_order12): SVD-Prony 分解返回频率、衰减、幅值、相位正频率部分 x np.asarray(x, dtypefloat) x x - x.mean() # 去直流避免零极点淹没低阶极点 N len(x) L N // 2 H np.array([x[i:i L] for i in range(N - L)]) U, S, Vh np.linalg.svd(H, full_matricesFalse) energy np.cumsum(S ** 2) / np.sum(S ** 2) p int(np.searchsorted(energy, 0.999)) 1 p max(4, min(p, max_order)) # 按奇异值能量定阶再设硬上限 V Vh[:p].T # (L, p) A np.linalg.lstsq(V[:-1], V[1:], rcondNone)[0] # 线性预测矩阵 poles np.linalg.eigvals(A) dt 1.0 / fs freqs np.angle(poles) / (2 * np.pi * dt) alphas np.log(np.abs(poles)) / dt keep (freqs 0) (np.abs(poles) 1.02) poles, freqs, alphas poles[keep], freqs[keep], alphas[keep] Z poles[None, :] ** np.arange(N)[:, None] coef, *_ np.linalg.lstsq(Z, x, rcondNone) amp, phase np.abs(coef), np.angle(coef) order np.argsort(freqs) return freqs[order], alphas[order], amp[order], phase[order]逻辑说明Hankel 矩阵的 SVD 用来定阶右奇异向量 V 构成信号子空间lstsq(V[:-1], V[1:])求相邻行之间的线性预测矩阵特征值就是 Prony 极点极点取对数得到频率和衰减再用最小二乘拟合每个极点的幅值和相位。参数说明max_order不建议超过 15阶数越高越容易过拟合L N//2是常用折中数据窗短时改到N//3|poles| 1.02的容差留给数值误差真正物理极点应当在单位圆上或圆内。3.3 DCMP 与 APM 的阈值选择与切换逻辑把极点分到直流、工频、主频三段之后论文定义了两个指标。DCMP 是直流分量幅值与直流主频幅值之比低合闸角时故障线路的衰减直流更重APM 是工频幅值归一化比与相位偏离平均相位程度的乘积高合闸角时故障线路的幅相都更突出。def dcmp_apm_select(signals, fs, dcmp_th0.7): n_lines signals.shape[0] dc_amp, fund_amp, fund_ph, main_amp [], [], [], [] for sig in signals: f, a, amp, ph prony_decompose(sig, fs) dc_amp.append(amp[np.argmin(np.abs(f))]) idx_f np.argmin(np.abs(f - 50)) fund_amp.append(amp[idx_f]) fund_ph.append(ph[idx_f]) others [i for i in range(len(f)) if abs(f[i] - 50) 5 and f[i] 5] main_amp.append(amp[others[np.argmax(amp[others])]] if others else amp[idx_f]) dcmp np.array(dc_amp) / (np.array(dc_amp) np.array(main_amp) 1e-12) ref np.mean(fund_ph) phase_dev np.abs(np.angle(np.exp(1j * (np.array(fund_ph) - ref)))) apm (np.array(fund_amp) / (np.max(fund_amp) 1e-12)) * (phase_dev / np.pi) if dcmp.max() dcmp_th: return int(np.argmax(dcmp)), DCMP, dcmp, apm return int(np.argmax(apm)), APM, dcmp, apm逻辑说明先按 Prony 结果取三个分量再分别计算 DCMP 和 APMdcmp.max() 0.7时走 DCMP 判据否则走 APM。参数说明dcmp_th0.7是论文仿真里的经验值现场要按消弧线圈补偿度和数据窗长度重整建议先跑 50~100 组录波做分布统计APM 里的相位参考用所有线路工频相位的平均如果有线路 CT 极性接反参考值会被污染。指标适用场景典型阈值主要失效模式DCMP低故障合闸角、衰减直流明显0.7数据窗太短直流没衰减完APM高故障合闸角、主频幅相差异大无固定阈取最大噪声抬高标准差相位参考被污染定阶 p两者共用4~12阶数过高出现虚假极点4. 消弧线圈并联中值电阻投切前后工频增量系数 K 选线4.1 等效电路与 K 值判据的成立条件消弧线圈并联中值电阻的系统正常运行时电阻短接或断开接地故障后短时投入电阻给故障线路叠加一个阻性电流。健全线路的零序电流主要由对地电容决定电阻投入前后变化有限故障线路的零序电流里多出一条经电阻回到中性点的通路工频分量明显增大。K 投切后工频幅值 / 投切前工频幅值故障线路 K 最大。成立条件投切前后零序电压基本不变CT 不饱和采样同步。消弧线圈的补偿状态要清楚如果是随动调谐投切前工频残流已经很小K 的分母太小会导致比值不稳定。常见做法是要求投切前工频幅值大于最小可辨识值比如二次侧 0.5 A否则放弃该次选线。4.2 投切时序与采样窗口的工程参数参数常见值说明采样率10 kHz仿真/ 4~12.8 kHz装置工频相量提取建议整周期 DFT投切前窗5 个周波避开故障初始暂态投切后等待2 个周波等电阻投入的暂态衰减投切后窗5 个周波与前置窗长度一致电阻值10~20 Ω 量级按系统电容电流和中性点电压整定最小工频电流0.5 A二次侧低于此值报数据无效投切时刻最好由装置同一时钟控制软件里根据开关量变位打时间戳避免自己猜。如果装置只能给出一路总启动信号投切前后窗的划分就要靠零序电压突变点做二次对齐。4.3 代码实现与抗干扰处理import numpy as np def incremental_k(before, after, fs10e3, f050.0): 投切前后工频增量系数before/after 形状 (n_lines, n_samples) def phasor(x): n int(len(x) // (fs / f0) * (fs / f0)) # 整周期窗 x x[:n] - np.mean(x[:n]) # 去均值压衰减直流偏置 k np.arange(n) ref np.exp(-1j * 2 * np.pi * f0 * k / fs) return 2.0 / n * np.sum(x * ref) k_values [] for i in range(before.shape[0]): i1 abs(phasor(before[i])) i2 abs(phasor(after[i])) # 分母太小说明投切前残流不可信置零并让上层报无效 k_values.append(i2 / i1 if i1 0.5 else 0.0) return np.array(k_values)逻辑说明用整周期 DFT 代替 FFT 求最大谱线避免非整周期泄漏去均值消掉衰减直流对工频相量的偏置分母加最小可信值判断。参数说明fs与装置采样率一致before窗和after窗长度相同0.5是二次侧门槛换成一次侧要乘 CT 变比。多线路 K 值排序后如果第二名的 K 超过第一名的 80%建议直接报“选线冲突”让运维结合暂态录波再判。注意K 值法对开关量时间戳很敏感投切前后窗如果错了一个周波健全线路的 K 也会被抬起来。5. 三方法横向验证与从 Matlab 录波到在线判线的迁移5.1 统一仿真测试框架与失败样本归因把三种方法放在同一组仿真信号上跑比单看论文里的曲线更有意义。常见做法是固定线路模型扫三个变量过渡电阻0.1 kΩ 到 5 kΩ、故障合闸角0° 到 90°、消弧线圈补偿度-10% 到 10%。每组变量下重复 50 次加 30 dB 噪声记录命中率。def sweep_test(make_signals, selectors, r_list, angle_list, trials50): make_signals 返回 (信号, 真实故障线路号)selectors 是方法名到函数的映射 result {} for name, func in selectors.items(): hit, total 0, 0 for r in r_list: for ang in angle_list: for _ in range(trials): sig, true_idx make_signals(r, ang) pred func(sig) hit int(pred true_idx) total 1 result[name] hit / total return result失败样本要按原因分类看相角差法失败大多出现在特征频率能量低于噪声底的情况Prony 法失败多出现在数据窗短于一个工频周期、或阶数给到 15 以上K 值法失败几乎都和投切前后窗不同步、CT 极性接反有关。5.2 从 Matlab 录波到 Python 在线判线的迁移技巧论文仿真和现场装置之间差在数据来源。Matlab 里常用load读.matPython 里对应scipy.io.loadmat要注意.mat里可能存的是single读进来先转float64否则lstsq容易报秩亏。通道顺序、采样率、CT 变比要单独记元数据不要靠文件名猜。滤波器群延迟必须补偿。前面用的 129 阶 FIR 在 10 kHz 采样下群延迟约 6.45 ms对应工频约 116°不补偿的话相角差判据直接失效。补偿方式算相位时把参考相位一起过同一个滤波器或者用scipy.signal.filtfilt做零相位滤波但filtfilt会改变暂态段K 值法不要用。如果装置只有 4 kHz 采样129 阶带通在 150~250 Hz 上过渡带会压到目标频带常见做法是把阶数降到 64 并接受更宽的过渡带或者改用二阶节级联的多相滤波。本文还有配套的精品资源点击获取 配电网小电流接地选线：零序导纳相角差、Prony与K值法复现
简介这份资源围绕小电流接地系统单相接地故障选线这一电力系统保护难题展开面向从事继电保护与控制研究的研究人员、工程师及研究生。针对传统选线方法在参数变化、高阻故障与不同合闸角下准确率偏低的问题资源系统复现了论文提出的三种新原理基于工频与特征频率零序导纳相角差判据、综合低频与主频分量的改进Prony分解方法以及消弧线圈并联中值电阻系统电阻投切前后工频增量系数法并给出配套代码与逐步解释。包内为1个pdf文件约910KB内含论文复现分析、方法原理推导、Matlab仿真流程与等效电路建模说明可对照代码理解特征频率搜索、DCMP与APM指标构建及故障线路判定逻辑。已有71人学习适合希望掌握选线判据实现与仿真验证、为装置选型和技术方案制定提供参考的读者。1. 小电流接地选线为什么在参数变化下失灵从三种新判据说起配电网里中性点不接地或经消弧线圈接地的系统单相接地时故障电流常常只有几安到几十安现场叫小电流接地系统。传统选线靠零序电流比幅比相最头疼的是过渡电阻一高、故障合闸角一小、消弧线圈补偿度一调幅值和相位跟着漂选线装置就开始乱报。这个复现包围绕论文里的三条新判据展开工频与特征频率零序导纳相角差、改进 Prony 分解后的 DCMP 与 APM 双判据以及消弧线圈并联中值电阻投切前后工频增量系数 K。代码以 Python 仿真为主配套 Matlab 风格的信号生成适合做保护算法、选线装置整定和论文复现的人拆着看。2. 零序导纳相角差判据特征频率搜索与相角差整定2.1 判据的物理依据与适用范围小电流接地系统的零序网络里健全线路的零序导纳主要来自对地电容工频下阻抗角接近 -90°故障线路的零序电流要流过接地点过渡电阻和消弧线圈高频段还叠加线路电感与对地电容构成的谐振分量。论文抓住的点是同一条线路在工频和特征频率下分别算零序导纳相角故障线路两者的角差接近 180°健全线路接近 0°。这条判据不依赖零序电流的绝对幅值对高阻故障比传统比幅法稳。它的边界也很明确特征频率必须真的存在。特征频率一般在 150~250 Hz 范围具体落点跟线路长度、电缆比例、消弧线圈补偿度有关如果系统里电缆很短、过渡电阻又特别大这段高频能量可能被 CT 和滤波器吃掉。另外各条线路 CT 的变比、极性必须一致否则相角差本身没有可比性。观测量健全线路典型值故障线路典型值整定注意工频零序导纳相角接近 -90°容性受消弧线圈影响上移取稳定窗平均特征频率零序导纳相角与工频接近与工频相差接近 180°特征频率要能量占优相角差绝对值40°120°门槛可先取 120°2.2 S 变换特征频率搜索的代码实现与坑原复现里的 S 变换写法window np.exp(-2*np.pi**2*f**2*t**2)直接拿t np.arange(N)/fs构造高斯窗等价于把窗中心放在 t0结果是对信号起始段加权最重对稳态段几乎没有分辨力。S 变换的窗中心应该随每个采样点移动工程上要快速复现常见做法是用带通滤波加滑动高斯窗做近似先找能量最大的候选频点。import numpy as np from scipy.signal import firwin, lfilter, hilbert def search_characteristic_freq(signals, fs, f_min150, f_max250, step5): 在候选频段内搜索总能量最大的频率signals 形状 (n_lines, n_samples) best_f, best_e None, -np.inf for f in np.arange(f_min, f_max step, step): # 带通取候选频点附近 ±5 Hz避免全带噪声参与能量统计 b firwin(129, [f - 5, f 5], pass_zeroFalse, fsfs) total 0.0 for sig in signals: y lfilter(b, 1, sig) env np.abs(hilbert(y)) # 希尔伯特包络近似该频点瞬时能量 total np.sum(env ** 2) if total best_e: best_f, best_e f, total return best_f, best_e逻辑说明每个候选频率做一次带通用希尔伯特包络平方累积能量哪个频点总能量最大就选为特征频率。参数说明fs是采样率现场录波常见 4 kHz、10 kHz、12.8 kHzf_min/f_max按系统电容电流水平先粗估10 kV 配网一般在 150~250 Hz35 kV 或长电缆系统可能要往上挪到 300 Hz滤波器阶数 129 对应 10 kHz 采样约 12.9 ms 窗长低于 2 个工频周期够用但会带来约 6.5 ms 群延迟后续算相角时必须做延迟补偿。提示能量搜索不要用原始信号的 FFT 峰值直接代替FFT 峰值容易落在工频泄漏的旁瓣上。2.3 工频与特征频率相角差的提取与判据整定知道特征频率之后接下来对每条线路分别提取工频和特征频率相量。稳妥做法是整周期 DFT 或窄带滤波加希尔伯特。def band_phase(sig, fs, f0, bw5): 提取 f0 附近分量的解缠相位 b firwin(129, [f0 - bw, f0 bw], pass_zeroFalse, fsfs) y lfilter(b, 1, sig) z hilbert(y) return np.unwrap(np.angle(z)) def phase_diff_select(signals, fs, f_char, f050): n_lines signals.shape[0] diff np.zeros(n_lines) for i, sig in enumerate(signals): th_f band_phase(sig, fs, f0) th_c band_phase(sig, fs, f_char) k len(sig) // 3 # 跳过滤波器启动段 d np.angle(np.exp(1j * (th_f[k:] - th_c[k:]))).mean() diff[i] abs(d) return diff, int(np.argmax(diff))逻辑说明band_phase先带通再希尔伯特np.angle(exp(1j*Δ))把相位差折到 (-π, π]避免直接相减出现 2π 跳变。参数说明bw5是半带宽频率成分密集时可缩到 2~3 Hz但滤波器阶数要跟着加k用来跳过启动段。判据整定建议最大角差超过 120° 才输出选线结果最大与次大差值小于 20° 时报“选线不确定”不要硬指一条。参数建议值调整方向特征频率搜索步长5 Hz频带拥挤时降到 2 Hz带通半带宽2~5 Hz噪声大取大分辨率要求高取小相角差门槛120°现场误报多就抬到 140°最大/次大间隔20°小于该值降级为不确定3. 改进 Prony 分解与 DCMP/APM 双判据实现3.1 为什么用 Prony 而不是直接 FFTFFT 受窗长限制10 kHz 采样、0.1 s 窗才 10 Hz 分辨率低频段和衰减直流分不开Prony 把信号建模成若干复指数之和能给出频率、阻尼、幅值和相位短窗更合适。小电流接地故障的零序电流里同时有工频残流、低频衰减分量和高频分量低合闸角时低频分量比高频更突出论文用 DCMP 抓这个特征高合闸角时主频分量的幅相差异更明显用 APM。这两种场景如果只用 FFT 的幅值比就容易被频率泄漏拉平。Prony 的代价是对噪声和模型阶数敏感阶数给大了会冒出虚假极点频率和阻尼都是假的。工程上不会直接套原始 Prony常见做法是 SVD 降阶加总体最小二乘或者先做带通再分解。3.2 Hankel 矩阵与极点求解的可运行写法原复现代码里V1 V[:, :n]、V2 V[:, 1:n1]、A pinv(V1[:-1]) V2[:-1]这一段下标逻辑是拧的V是右奇异向量构成的矩阵按那个写法求出的 A 既不是方阵也没有明确的极点含义。更接近标准 SVD-Prony 的写法如下。import numpy as np def prony_decompose(x, fs, max_order12): SVD-Prony 分解返回频率、衰减、幅值、相位正频率部分 x np.asarray(x, dtypefloat) x x - x.mean() # 去直流避免零极点淹没低阶极点 N len(x) L N // 2 H np.array([x[i:i L] for i in range(N - L)]) U, S, Vh np.linalg.svd(H, full_matricesFalse) energy np.cumsum(S ** 2) / np.sum(S ** 2) p int(np.searchsorted(energy, 0.999)) 1 p max(4, min(p, max_order)) # 按奇异值能量定阶再设硬上限 V Vh[:p].T # (L, p) A np.linalg.lstsq(V[:-1], V[1:], rcondNone)[0] # 线性预测矩阵 poles np.linalg.eigvals(A) dt 1.0 / fs freqs np.angle(poles) / (2 * np.pi * dt) alphas np.log(np.abs(poles)) / dt keep (freqs 0) (np.abs(poles) 1.02) poles, freqs, alphas poles[keep], freqs[keep], alphas[keep] Z poles[None, :] ** np.arange(N)[:, None] coef, *_ np.linalg.lstsq(Z, x, rcondNone) amp, phase np.abs(coef), np.angle(coef) order np.argsort(freqs) return freqs[order], alphas[order], amp[order], phase[order]逻辑说明Hankel 矩阵的 SVD 用来定阶右奇异向量 V 构成信号子空间lstsq(V[:-1], V[1:])求相邻行之间的线性预测矩阵特征值就是 Prony 极点极点取对数得到频率和衰减再用最小二乘拟合每个极点的幅值和相位。参数说明max_order不建议超过 15阶数越高越容易过拟合L N//2是常用折中数据窗短时改到N//3|poles| 1.02的容差留给数值误差真正物理极点应当在单位圆上或圆内。3.3 DCMP 与 APM 的阈值选择与切换逻辑把极点分到直流、工频、主频三段之后论文定义了两个指标。DCMP 是直流分量幅值与直流主频幅值之比低合闸角时故障线路的衰减直流更重APM 是工频幅值归一化比与相位偏离平均相位程度的乘积高合闸角时故障线路的幅相都更突出。def dcmp_apm_select(signals, fs, dcmp_th0.7): n_lines signals.shape[0] dc_amp, fund_amp, fund_ph, main_amp [], [], [], [] for sig in signals: f, a, amp, ph prony_decompose(sig, fs) dc_amp.append(amp[np.argmin(np.abs(f))]) idx_f np.argmin(np.abs(f - 50)) fund_amp.append(amp[idx_f]) fund_ph.append(ph[idx_f]) others [i for i in range(len(f)) if abs(f[i] - 50) 5 and f[i] 5] main_amp.append(amp[others[np.argmax(amp[others])]] if others else amp[idx_f]) dcmp np.array(dc_amp) / (np.array(dc_amp) np.array(main_amp) 1e-12) ref np.mean(fund_ph) phase_dev np.abs(np.angle(np.exp(1j * (np.array(fund_ph) - ref)))) apm (np.array(fund_amp) / (np.max(fund_amp) 1e-12)) * (phase_dev / np.pi) if dcmp.max() dcmp_th: return int(np.argmax(dcmp)), DCMP, dcmp, apm return int(np.argmax(apm)), APM, dcmp, apm逻辑说明先按 Prony 结果取三个分量再分别计算 DCMP 和 APMdcmp.max() 0.7时走 DCMP 判据否则走 APM。参数说明dcmp_th0.7是论文仿真里的经验值现场要按消弧线圈补偿度和数据窗长度重整建议先跑 50~100 组录波做分布统计APM 里的相位参考用所有线路工频相位的平均如果有线路 CT 极性接反参考值会被污染。指标适用场景典型阈值主要失效模式DCMP低故障合闸角、衰减直流明显0.7数据窗太短直流没衰减完APM高故障合闸角、主频幅相差异大无固定阈取最大噪声抬高标准差相位参考被污染定阶 p两者共用4~12阶数过高出现虚假极点4. 消弧线圈并联中值电阻投切前后工频增量系数 K 选线4.1 等效电路与 K 值判据的成立条件消弧线圈并联中值电阻的系统正常运行时电阻短接或断开接地故障后短时投入电阻给故障线路叠加一个阻性电流。健全线路的零序电流主要由对地电容决定电阻投入前后变化有限故障线路的零序电流里多出一条经电阻回到中性点的通路工频分量明显增大。K 投切后工频幅值 / 投切前工频幅值故障线路 K 最大。成立条件投切前后零序电压基本不变CT 不饱和采样同步。消弧线圈的补偿状态要清楚如果是随动调谐投切前工频残流已经很小K 的分母太小会导致比值不稳定。常见做法是要求投切前工频幅值大于最小可辨识值比如二次侧 0.5 A否则放弃该次选线。4.2 投切时序与采样窗口的工程参数参数常见值说明采样率10 kHz仿真/ 4~12.8 kHz装置工频相量提取建议整周期 DFT投切前窗5 个周波避开故障初始暂态投切后等待2 个周波等电阻投入的暂态衰减投切后窗5 个周波与前置窗长度一致电阻值10~20 Ω 量级按系统电容电流和中性点电压整定最小工频电流0.5 A二次侧低于此值报数据无效投切时刻最好由装置同一时钟控制软件里根据开关量变位打时间戳避免自己猜。如果装置只能给出一路总启动信号投切前后窗的划分就要靠零序电压突变点做二次对齐。4.3 代码实现与抗干扰处理import numpy as np def incremental_k(before, after, fs10e3, f050.0): 投切前后工频增量系数before/after 形状 (n_lines, n_samples) def phasor(x): n int(len(x) // (fs / f0) * (fs / f0)) # 整周期窗 x x[:n] - np.mean(x[:n]) # 去均值压衰减直流偏置 k np.arange(n) ref np.exp(-1j * 2 * np.pi * f0 * k / fs) return 2.0 / n * np.sum(x * ref) k_values [] for i in range(before.shape[0]): i1 abs(phasor(before[i])) i2 abs(phasor(after[i])) # 分母太小说明投切前残流不可信置零并让上层报无效 k_values.append(i2 / i1 if i1 0.5 else 0.0) return np.array(k_values)逻辑说明用整周期 DFT 代替 FFT 求最大谱线避免非整周期泄漏去均值消掉衰减直流对工频相量的偏置分母加最小可信值判断。参数说明fs与装置采样率一致before窗和after窗长度相同0.5是二次侧门槛换成一次侧要乘 CT 变比。多线路 K 值排序后如果第二名的 K 超过第一名的 80%建议直接报“选线冲突”让运维结合暂态录波再判。注意K 值法对开关量时间戳很敏感投切前后窗如果错了一个周波健全线路的 K 也会被抬起来。5. 三方法横向验证与从 Matlab 录波到在线判线的迁移5.1 统一仿真测试框架与失败样本归因把三种方法放在同一组仿真信号上跑比单看论文里的曲线更有意义。常见做法是固定线路模型扫三个变量过渡电阻0.1 kΩ 到 5 kΩ、故障合闸角0° 到 90°、消弧线圈补偿度-10% 到 10%。每组变量下重复 50 次加 30 dB 噪声记录命中率。def sweep_test(make_signals, selectors, r_list, angle_list, trials50): make_signals 返回 (信号, 真实故障线路号)selectors 是方法名到函数的映射 result {} for name, func in selectors.items(): hit, total 0, 0 for r in r_list: for ang in angle_list: for _ in range(trials): sig, true_idx make_signals(r, ang) pred func(sig) hit int(pred true_idx) total 1 result[name] hit / total return result失败样本要按原因分类看相角差法失败大多出现在特征频率能量低于噪声底的情况Prony 法失败多出现在数据窗短于一个工频周期、或阶数给到 15 以上K 值法失败几乎都和投切前后窗不同步、CT 极性接反有关。5.2 从 Matlab 录波到 Python 在线判线的迁移技巧论文仿真和现场装置之间差在数据来源。Matlab 里常用load读.matPython 里对应scipy.io.loadmat要注意.mat里可能存的是single读进来先转float64否则lstsq容易报秩亏。通道顺序、采样率、CT 变比要单独记元数据不要靠文件名猜。滤波器群延迟必须补偿。前面用的 129 阶 FIR 在 10 kHz 采样下群延迟约 6.45 ms对应工频约 116°不补偿的话相角差判据直接失效。补偿方式算相位时把参考相位一起过同一个滤波器或者用scipy.signal.filtfilt做零相位滤波但filtfilt会改变暂态段K 值法不要用。如果装置只有 4 kHz 采样129 阶带通在 150~250 Hz 上过渡带会压到目标频带常见做法是把阶数降到 64 并接受更宽的过渡带或者改用二阶节级联的多相滤波。本文还有配套的精品资源点击获取