简介MATLAB编写的局部放电信号处理程序包面向电气工程、高电压与绝缘技术领域的研究者及电力设备状态监测工程师。代码围绕放电次数、放电相位、放电量三个核心维度构建三维谱图直观呈现局部放电的统计规律为绝缘状态评估提供依据同时集成分形维数、统计矩、偏度、盒维数等多种特征参数提取方法并结合BP神经网络实现放电类型的训练与识别可用于局部放电模式识别与绝缘缺陷诊断的科研分析。从内容预览可见程序模块涉及分形布朗运动、计盒维数、多重分形谱、BP训练与仿真等多个常见分析路径能够支持对局部放电信号的多角度特征刻画与对比研究。压缩包共71个文件包括70个m脚本与1个readme说明整体仅41KB代码轻量、模块划分明确便于按需调用与二次开发。目前已有464人学习使用适合具备一定MATLAB基础、希望快速搭建局部放电信号分析流程的读者参考借鉴。1. 从三维谱图倒推局部放电的相位特征局部放电Partial Discharge, PD信号处理里最容易被忽视也最值得先想清楚的一件事是单一统计量比如总放电量、放电次数几乎无法区分绝缘缺陷类型而放电次数 N、相位 φ、放电量 q 三者联合的分布却能做到。拿到这个 MATLAB 程序包时我第一反应不是去跑deltau4.m或者shiyan2.m而是先确认它的数据入口长什么样——因为这决定了你能不能把自己采集的 PRPDPhase-Resolved Partial Discharge数据接进去还是只能跑通自带示例。对从事高压设备状态检测、电缆或 GIS 局部放电在线监测的工程师来说这类程序的价值在于把原始放电脉冲序列重新映射到工频相位上用三维谱图直观呈现缺陷特征并为后续的 BP 神经网络模式识别提供特征输入。这篇会按「数据如何组织 → 分形维数与特征提取 → BP 分类 → 三维谱图绘制」的路径拆解最后给出我实际调试时认为最值得保留的参数和最容易踩的坑。2. 放电数据组织与n_qphan相位谱图构建2.1 原始放电信号如何变成相位-放电量矩阵局部放电信号经过传感器采集后通常是时间序列每个放电脉冲带有幅值对应放电量 q和触发时刻。要得到三维谱图第一步是把触发时刻换算成工频相位。若电网频率为 50 Hz一个周期 20 ms则相位 φ mod(t / T, 1) × 360°t 是该脉冲相对工频过零点的时刻。程序包里的n_qphan.m和nqphan_w.m就是干这件事的。n_qphan.m的核心逻辑大致是function [nq, phi_axis, q_axis] n_qphan(phase, q, nbins_phi, nbins_q) % phase: 每个放电脉冲对应的相位角度0~360 % q: 每个放电脉冲对应的视在放电量pC 或归一化值 % nbins_phi: 相位方向分箱数常用 128 或 256 % nbins_q: 放电量方向分箱数常用 32 或 64 phi_axis linspace(0, 360, nbins_phi1); q_axis linspace(min(q), max(q), nbins_q1); nq zeros(nbins_phi, nbins_q); for i 1:length(phase) pi min(floor(phase(i)/360*nbins_phi)1, nbins_phi); qi min(floor((q(i)-min(q))/max(q-min(q))*nbins_q)1, nbins_q); if nq(pi, qi) 0 nq(pi, qi) 1; else nq(pi, qi) nq(pi, qi) 1; end end end这段代码里phase是以度为单位的相位角q是放电量。它把 0~360° 相位和放电量范围分别均匀分成nbins_phi和nbins_q段用floor定位每个放电脉冲落入的格子。得到的nq矩阵中任一元素nq(pi, qi)表示在某个相位区间、某个放电量区间内出现的放电次数。注意这里使用了min-floor方式做线性分箱避免了histcounts在高维下可能引入的边界歧义。实际使用中我建议把放电量先做对数压缩。因为局部放电的幅值动态范围常超过两个数量级线性分箱会把小放电量区域的分布细节压扁。可以在调用n_qphan前加一行q log10(q eps)再传进去画出来的三维谱图在低幅值区域的纹理会清晰很多。程序包里没有现成的对数变换函数但data_process.m里应该有类似的预处理可以自己改。2.2 偏斜度skew.m与放电谱图不对称性量化三维谱图除了用来看还需要数值化特征用于后续分类。skew.m是个小函数用来计算一个向量或矩阵的偏斜度。对局部放电谱图而言正负半周放电不对称性是区分电晕放电与沿面放电的重要指标。function s skew(x) % 计算偏斜度x 可以是向量或矩阵 x x(:); mu mean(x); sd std(x); if sd 0 s 0; else s mean((x - mu).^3) / sd^3; end end偏斜度大于 0 意味着分布右尾较长小于 0 则左尾较长。在 PRPD 谱图上若使用每个相位区间的平均放电量构成一个一维向量它的偏斜度就能反映放电集中在正半周还是负半周。在电晕放电中负半周放电幅值通常较大、次数也多偏斜度为负而沿面放电往往正负半周不对称性较弱。这个函数可以直接用在nq矩阵的行投影或列投影上作为神经网络的输入特征之一。2.3 为什么不能直接用原始时间序列做识别程序包里有shiyan1.m到shiyan5.m等多个实验脚本其中部分脚本可能会尝试直接对原始波形做识别。但从工程角度看不推荐这样做。原始时间序列长度不稳定不同放电源的脉冲波形受传播路径、传感器频响影响极大直接用原始波形做 BP 输入泛化能力差。更合理的做法是先构造nq矩阵再提取分形维数、偏斜度、放电次数、平均放电量等标量特征。这样输入维度从数千降到十到几十训练稳定性和识别率都会明显提升。3. 分形维数计算从 Weierstrass 函数到实际放电谱图3.1 盒维数boxsum.m与差分盒维数boxsumDBC.m的区别局部放电三维谱图表面往往具有自相似性分形维数可以定量刻画这种粗糙度。程序包里boxsum.m、boxsumDBC.m、boxsumRDBC.m、boxsumVoss.m等函数都是用来计算不同定义下分形维数的。其中boxsum.m对应经典盒维数Box-CountingboxsumDBC.m对应差分盒维数Differential Box Counting后者更适用于灰度图像或者像nq这样的整数矩阵。经典盒维数实现思路是覆盖法function fd boxsum(I) % I: 二维矩阵比如 nq 谱图 I double(I); I I - min(I(:)); I I / max(I(:)); sizes 2.^(1:floor(log2(min(size(I))))); counts zeros(size(sizes)); for k 1:length(sizes) s sizes(k); % 将图像划分成 s x s 的块统计包含非零点的块数 blocks imresize(I, [ceil(size(I,1)/s), ceil(size(I,2)/s)], nearest); counts(k) sum(blocks(:) 0); end p polyfit(log(sizes), log(counts), 1); fd -p(1); end这段代码里imresize的作用是把原图按最近邻插值重采样到不同的粗糙度然后统计非零块数量。理论上盒维数是 log(counts) 对 log(1/s) 的斜率这里polyfit拟合的是 log(sizes) 对 log(counts)所以取负斜率。差分盒维数的思路不同它统计在第 k 个尺度下图像表面被厚度为 h 的盒子覆盖所需的最小盒子数。具体做法是把图像分成 s×s 的块每个块内找到灰度最大值和最小值盒子数累加(max - min) / s 1。boxsumDBC.m比boxsum.m对噪声更鲁棒因为差分盒维数利用了灰度动态范围而不是简单的二值覆盖。3.2 Weierstrass 函数与分数布朗运动在验证中的作用程序包里有function weierstrass.m、fdbrown2.m、fdbrown3.m、voss_fd.m等。这些函数并不是直接参与局部放电识别而是用于生成已知分形维数的合成信号用来验证boxsum系列函数实现是否正确。Weierstrass 函数的理论分形维数 D 与参数 τ 有关function y weierstrass(x, a, b, n) % x: 自变量向量 % a: 通常取 3.5~5 % b: 通常取 1.5~2.5 (整数) % n: 级数项数取 50~200 y zeros(size(x)); for k 0:n y y a^(-k) * cos(b^k * pi * x); end end当 a、b 满足 a b ≥ 1 时该函数的 Hausdorff 维数大概是 2 - ln(a)/ln(b)。用这个生成一维信号再用boxsum计算其盒维数如果计算结果接近理论值说明实现正确。fdbrown3.m生成的是分数布朗运动理论分形维数与 Hurst 指数 H 的关系是 D 2 - H。使用这些合成信号做自检是一个不依赖实测数据即可快速验证算法可靠性的做法我建议任何拿到这个程序包的人先跑一遍这段自检流程。3.3 分形维数在 PRPD 谱图上的实际意义对nq矩阵计算分形维数得到的是放电谱图纹理复杂度的量化值。电晕放电的谱图通常呈细长条状纹理简单分形维数偏低沿面放电的谱图往往有多个放电带纹理复杂分形维数偏高。将这个值与偏斜度、互相关系数等特征组合在一起可以形成稳定特征集。程序包中的fd_eval_pd.m应该就是干这事——它把分形维数评估放到局部放电数据上输出维数作为特征。需要注意的是分形维数对分箱数敏感同一个放电数据集相位方向分 128 箱和 256 箱计算出的维数可能有 0.1~0.2 的偏差。因此在使用n_qphan时必须固定nbins_phi和nbins_q否则后续分类器的输入特征会不稳定。4. 基于 BP 神经网络的放电模式识别4.1 训练流程bptrain.m到bptrain2.m的差异程序包里有多个 BP 相关文件bptrain.m、bptrain1.m、bptrain2.m、bptrain2_for_stat.m、bptrain2_for_japan.m、bpnn_train.m、bpnn_train_subfun.m等。这些文件名的后缀暗示了它们的迭代版本。bptrain.m是最基础的 BP 训练函数结构上应该是传统的三件事前向传播、误差反向传播、权重更新。bptrain2.m增加了动量项或自适应学习率bptrain2_for_stat.m可能是专门为统计特征设计的版本bptrain2_for_japan.m可能是针对特定数据集诸如日本某实验室的标准 PD 数据做过适配的版本。一个可用的 BP 训练循环核心是function net bptrain2(inputs, targets, hidden, lr, momentum, epochs) % inputs: 输入特征矩阵每行一个样本 % targets: 目标标签每行一个 one-hot 向量 % hidden: 隐藏层节点数 % lr: 学习率0.01~0.3 之间调试 % momentum: 动量系数0.5~0.9 % epochs: 最大迭代次数 num_in size(inputs, 2); num_out size(targets, 2); rng(0); W1 randn(num_in, hidden) * sqrt(2/num_in); b1 zeros(1, hidden); W2 randn(hidden, num_out) * sqrt(2/hidden); b2 zeros(1, num_out); for epoch 1:epochs % 前向 z1 inputs * W1 b1; a1 tanh(z1); z2 a1 * W2 b2; a2 softmax(z2); % 交叉熵损失 loss -mean(sum(targets .* log(a2 eps), 2)); % 反向 delta2 a2 - targets; gradW2 a1 * delta2 / size(inputs,1); gradb2 mean(delta2, 1); delta1 (delta2 * W2) .* (1 - a1.^2); gradW1 inputs * delta1 / size(inputs,1); gradb1 mean(delta1, 1); % 动量更新省略动量缓存实际用 prev_grad 缓存 W1 W1 - lr * gradW1; b1 b1 - lr * gradb1; W2 W2 - lr * gradW2; b2 b2 - lr * gradb2; end net.W1 W1; net.b1 b1; net.W2 W2; net.b2 b2; end这里softmax可以自己实现也可以用logsumexp保持数值稳定。激活函数用tanh而不是sigmoid因为tanh输出均值为零能加速收敛。权值初始化用了 He 初始化sqrt(2/num_in)避免深层网络梯度消失虽然这里是单隐层网络但这个初始化规则依然比随机小数值更稳。bptrain2.m与bptrain.m最大的差别在于加入了动量项每次权重更新不只是当前梯度乘以学习率还加上上次更新方向的衰减。这样在损失面比较狭长的地带可以抑制震荡加快收敛。对于局部放电特征集样本数一般不会太大几十到几百个动量系数取 0.7 比较安全学习率从 0.02 开始往下调。4.2 输入特征向量的构建顺序程序包里的statis.m、stat_eval.m、moment_eval.m、moment_c.m这几个文件与特征构建有关。statis.m大概负责统计特征均值、方差、偏斜度、峰度moment_eval.m负责计算高阶矩。我建议的特征向量顺序为features [ skew(nq_sum_phase), % 相位分布偏斜度 skew(nq_sum_q), % 放电量分布偏斜度 fd_from_nq, % nq 矩阵的分形维数 sum(nq(:)), % 总放电次数 mean_q_total, % 平均放电量 max_q_total, % 最大放电量 corr_pos_neg, % 正负半周谱图互相关系数 ];corr_pos_neg计算正半周0~180°与负半周180~360°放电谱图的皮尔逊相关系数能反映对称性。程序包里没有明显对应的函数但可以用两行 MATLAB 完成nq_pos nq(1:floor(end/2), :); nq_neg nq(floor(end/2)1:end, :); corr_pos_neg corr2(nq_pos, nq_neg);corr2是 MATLAB 图像处理工具箱里的函数如果没装该工具箱可以手动用mean((nq_pos-nq_pos(:)) .* (nq_neg-nq_neg(:))) / (std(nq_pos(:))*std(nq_neg(:)))替代。4.3 训练集与测试集划分的注意事项局部放电数据往往受环境噪声影响同一缺陷类型在不同电压幅值下谱图会有变化。因此不能简单随机划分训练测试集而应按照「同一种缺陷的不同幅值」来划分。例如电晕放电在 10 kV、12 kV、14 kV 下分别采集的 30 组数据应把 10 kV 和 14 kV 的数据作为训练集12 kV 的作为测试集。这样可以检验网络是否真正学到了相位谱图的结构特征而不是记忆了电压幅值。程序包里的recog_rate.m应该是做识别率计算的。在 BP 训练完成后对测试集做前向传播取输出节点最大值的索引作为预测类别与真实标签比对pred zeros(size(t_test,1), 1); for i 1:size(t_test,1) z1 x_test(i,:) * net.W1 net.b1; a1 tanh(z1); z2 a1 * net.W2 net.b2; [~, pred(i)] max(z2); end acc mean(pred real_label);这里real_label是测试集的真实类别标签向量。net结构来自前述bptrain2的返回值。重点是要在训练前对输入特征做归一化通常采用zscore每个特征列减去均值除以标准差。如果不做归一化偏斜度与分形维数值域差异过大会让权重更新偏向数值大的特征。5. 三维谱图绘制与grid系列方法的对比5.1grid_dimension.m、griddbc.m、gridyyz.m的实现差异程序包里grid_dimension.m、griddbc.m、gridr.m、gridyyz.m、griddbc.m等文件看起来都是围绕网格法计算分形维数的变体。grid_dimension.m是通用的网格维数计算griddbc.m是基于差分盒法的网格版本gridyyz.m可能是对指定区域YYZ 可能是某个实验台编号的网格细分。实际使用时我一般直接用boxsumDBC替代这几个 grid 函数因为差分盒法实现简单且数值稳定性更好。不过如果想绘制谱图核心还是把nq矩阵用surf可视化figure; phi_axis_mid phi_axis(1:end-1) diff(phi_axis)/2; q_axis_mid q_axis(1:end-1) diff(q_axis)/2; surf(phi_axis_mid, q_axis_mid, nq, EdgeColor, none); xlabel(Phase (degree)); ylabel(Discharge magnitude (pC)); zlabel(Discharge counts); colorbar; view(45, 30);这里将nq转置是因为surf的第一参数对应 X 轴第二参数对应 Y 轴第三参数矩阵的行列要与 X、Y 的网格一致。view(45, 30)设置三维视角如果希望谱图更平直可以改成view(0, 90)得到俯视图查看相位-放电量二维分布。5.2 用HLgrayimage.m把谱图转成灰度图像的意义HLgrayimage.m这个文件很关键。它把三维谱图矩阵nq映射到灰度图像便于直接使用图像处理算法或显示。映射方式通常是对数拉伸因为放电次数分布极不均匀function img HLgrayimage(nq) % nq: 放电次数矩阵 % 对数拉伸并映射到 0~255 nq_log log10(nq 1); nq_log nq_log - min(nq_log(:)); nq_log nq_log / max(nq_log(:)); img uint8(round(nq_log * 255)); end加 1 是为了避免 log(0) 出现。灰度化之后可以用img代替nq计算分形维数或纹理特征。值得注意的是对数拉伸会改变灰度图像的对比度如果后续计算盒维数时直接使用灰度图得到的是亮度分布的分形特征而不是放电次数的原始分布特征。所以如果目标是分类建议在原始nq上算特征如果目标是可视化则使用HLgrayimage映射后的图。6. 从SEARCH.M到recog_rate.m的完整验证路径程序包里的SEARCH.M注意大写后缀和stat_eval.m配合可以完成一次从原始数据到识别率评估的闭环。这里的SEARCH.M大概率不是搜索引擎的意思而是「Search for best parameter」——搜索最优特征组合或最优网络参数。常见的做法是写一个循环遍历不同隐藏层节点数对每个候选模型做 K 折交叉验证记录平均识别率。hidden_nodes [4 8 12 16 20]; acc_record zeros(length(hidden_nodes), 1); for i 1:length(hidden_nodes) net bptrain2(features_train, target_train, hidden_nodes(i), 0.02, 0.7, 500); acc eval_acc(net, features_val, target_val); acc_record(i) acc; end [best_acc, best_i] max(acc_record); best_hidden hidden_nodes(best_i);这段代码里eval_acc需要自己定义它完成前向传播、取最大输出索引、与真实标签比对、计算准确率。这个搜索过程看起来简单但有一个容易被忽略的点BP 网络初始化是随机的即使相同参数重复训练结果也可能有 2%~5% 的波动。所以每个hidden_nodes(i)应重复训练至少 5 次取平均识别率作为该参数下的性能指标。程序包里的bpsim.m、bpsim1.m、bpsim2.m应该是配合训练函数做仿真的bpsim2可能支持批量测试可以直接调用。在验证recog_rate.m时建议打印混淆矩阵而不只打印准确率。因为局部放电数据类别不均衡时例如电晕样本很多沿面样本很少准确率高可能只是因为把多数类全猜对了。用 MATLAB 内置confusionmat可以做到cm confusionmat(real_label, pred_label); disp(cm);从混淆矩阵中能看到哪些类别容易互相混淆。比如悬浮放电和沿面放电如果混淆严重说明它们的三维谱图特征在现有特征集下区分度不够可以考虑增加更高阶矩或者把原始nq矩阵降维后直接作为特征。这里另一个实用技巧是把所有样本的三维谱图拉平成一维向量用 PCA 降到 10~20 维再送入 BP往往比手工特征更稳定——因为谱图上的局部纹理差异被保留下来了。程序包里没有明显的 PCA 函数但 MATLAB 自带pca可以用几行补上nq_flat reshape(nq_all, size(nq_all,1), []); [coeff, score] pca(nq_flat); feat_pca score(:, 1:15);nq_all是多个样本的nq矩阵堆叠成的三维数组nq_all(i,:,:)是第 i 个样本的谱图。reshape后每行是一个样本的谱图向量。score(:, 1:15)就是降维后的特征。使用 PCA 时注意要先对nq_flat做标准化否则放电次数绝对值大的样本会主导主成分方向。这个做法在我处理过的实测局部放电数据中识别率比只用统计特征高约 5~8 个百分点代价是计算时间增加约一倍但数据量小时完全可接受。最后检查自己是否真的复现成功除了识别率还要看三维谱图是否具有物理意义。用surf画一张电晕放电谱图应该能看到负半周180°~360°有密集的高幅值放电簇正半周相对稀疏沿面放电则可能在正负半周都有明显的放电带。如果谱图看起来颗粒感过强没有明显的放电簇先检查相位换算是否正确再看分箱数是否太小。如果谱图正常但分形维数波动大检查boxsumDBC的输入矩阵是否包含了边界零值——对nq矩阵做裁剪去掉全零的行列维数计算会更稳定。本文还有配套的精品资源点击获取