基于变分模态分解和凹凸型阈值小波的电缆局部放电信号降噪方法

2022-10-17吴昊王东山

现代电力 2022年5期

吴昊，王东山

（1. 上海电力大学电子与信息工程学院, 上海市杨浦区 200090；2. 北京智芯微电子科技有限公司, 北京市昌平区 102200）

0 引言

局部放电(partial discharge，PD)是电力电缆绝缘老化、劣化的表现形式之一，同时也是绝缘进一步劣化导致绝缘失效的主要原因之一[1]。所以PD监测是电缆绝缘状态监测、绝缘故障定位的重要手段。电缆局部放电信号常用的在线监测方法有脉冲电流法、超声波法、特高频法[2]。而在实际监测中，获得的实测信号往往存在大量的噪声干扰，而随机脉冲型干扰[3]、周期性窄带干扰[4-5]、白噪声干扰[6]是实测信号中较为常见的噪声类型。这些噪声问题大大降低了电缆局部放电在线监测的准确性。本文主要研究如何有效抑制局部放电信号中的周期性窄带干扰、白噪声干扰。

PD信号是一种短时突变的非平稳信号[7]。针对PD监测中的噪声问题常用的降噪方法有小波分析、经验模态分解(empirical mode decomposition，EMD)、变分模态分解(variational mode decomposition, VMD)等。小波降噪效果和分解尺度、阈值等的选择密切相关，不好的选择会在降噪过程中产生偏差、畸变、特征信息损失等情况。EMD缺点是无法确定PD信号中干扰信号的频带中心与宽度，存在模态混叠[8]、端点效应等问题，所以整体的降噪效果一般。而VMD可以将噪声信号和PD信号分解到几个独立的频带内，再通过信号重构达到抑制噪声的目的[9]。这样就摆脱了EMD中存在的模态混叠问题，但 VMD 并不能一步到位实现对非平稳信号的分离[10-11]。使用VMD分解重构后的有效信号中仍然会存在噪声残留。针对VMD噪声残留问题文献[10-11]中均使用传统的小波阈值方法来进行进一步的降噪处理，而传统的小波阈值在降噪过程存在上文提到的小波降噪缺陷。

针对以上降噪方法中存在的问题，本文提出了一种基于VMD和凹凸型阈值的小波降噪相结合的PD信号降噪方法。首先信号经过VMD分解重构，再使用凹凸型阈值的小波降噪方法完成对信号的降噪。使用本文提出的方法对PD仿真信号和实测信号进行降噪处理，并与小波降噪、VMD和传统小波相结合的方法进行对比，验证本文提出的方法在降噪效果上的提升。

1 基于VMD和凹凸型阈值小波变换的降噪方法

1.1 VMD初步降噪

1.1.1 VMD分解原理

本文首先使用VMD对信号进行分解重构来达到信号初步降噪的目的。相较于通过循环迭代逐次分解本征模态（intrinsic mode function, IMF）的经验模态分解。VMD核心思想就是通过构建和求解变分问题来完成对原始信号f(t)的模态分解，得到K个以ωk为中心频率的模态uk。具体的变分模型表达式如下

式中： {uk}={u1,···,uK}，{ωk}={ω1,···,ωK}；K表示分解的IMF的总个数；表示通过 Hibert变换和卷积运算得到uk(t)的解析信号；e-jωkt表示将每个模态uk(t)的解析信号频谱平移到基带上去；最后通过计算平方L2范数来估计各模态的带宽。通过引入二次惩罚因子α和Lagrange乘子λ(t)转换为非约束的形式如公式（2）所示

利用交替方向乘子法对公式（2）进行求解，在频域对uk,ωk和λ进行更新迭代。迭代过程如式（3）（4）（5）所示：

1.1.2 峭度准则

峭度用来表示信号分布特性的归一化4阶数值统计量，它是一种无量纲参数，用来表示波形的尖峰度[12-13]。对于一个离散信号，峭度的表达式为

式中：µ表示信号均值；σ表示信号标准差；E(x-µ)4表示4阶数学期望。PD信号作为一种短时、陡峭的上升信号，含有PD信号的模态将有更高的峭度值。本文将使用峭度值来确定VMD分解的K值、参与重构的模态、惩罚因子 α值。

1.1.3 VMD参数调试

使用VMD进行分解过程中主要需要调试的内容有K值确定、重构模态的选择、α值的确定，调试主要内容如下。

1）K值的确定。

当K值分解数过大时会存在过度分解现象，所以本文的K值选择范围为[2,7]。在仿真实验阶段K值可以通过分解后的模态中心频率和PD信号的频率是否接近来选择合适的K值。而对于实测的PD信号，其信号的中心频率属于未知量，可以通过模态平均峭度的大小来选择合适的K值，平均峭度值最高的K值最合适。本文在仿真实验中将分别使用2种方法来确定K值，对于实测信号将使用计算平均峭度来确定K值。

2）参与重构的模态选择。

参与信号重构的模态选择也需要通过峭度准则来判断。不含有PD信号特征的模态峭度值在3左右，含有PD信号特征的模态峭度值将远大于3[12]。

3）惩罚因子 α的确定。

惩罚因子α决定了VMD分解中各模态的寻优频段范围，α值越大则频段范围越小，反之越小则频段范围越大。所以在K值固定的情况下随着α值的增大，模态的寻优频段范围将不断接近PD信号的频段范围，这个阶段平均峭度值和信噪比都将不断增加。而当α值增大到一定程度时，寻优频段范围将会小于PD信号的频段范围，此时PD信号特征的损失将不断增加，峭度值将会停止增加甚至开始减少。所以可以通过寻找模态平均峭度取得最大值时 α值作为VMD分解参数。

1.2 凹凸型阈值小波变换降噪

由于白噪声的带宽分布在全频段内，所以经过VMD分解重构后的信号仍有与PD信号处在相同或相近频带的噪声残留。所以这里需要小波降噪方法对重构后的信号进行进一步的降噪。

硬阈值函数式(7)和软阈值函数[14]式(8)的函数公式如下所示

软硬阈值函数在阈值的边界处会存在不连续性和偏差的问题，为了解决这些问题，本文引入一种全新的改进阈值函数[15]，使用凹度和指数函数的凸性来重构的凹凸型阈值函数为

式中：ws(j,k)表示各层的小波系数；j表示分解的尺度；k表示小波层数；λ为小波阈值；m（m>=1）取正整数为指数调整系数，这里取默认值2。3种阈值函数的对比如图1所示。

图1 阈值函数的对比Fig. 1 Comparison of threshold functions

本文引入自适应阈值[16]，公式如下

式中：N表示信号的长度；σ为信号标准差，标准差计算公式如公式(11)所示。自适应阈值大小和小波分解规模呈负相关，而不是传统的固定阈值。这样可以更有效地区分原始信号的噪声，提升降噪效果。

1.3 降噪方法的流程

对于实际PD信号的处理流程如下所示。

1）使用α= 2000、τ=0作为VMD分解的初始参数，K值的取值范围设为[2,7]。对输入的原始信号进行VMD分解，取模态平均峭度最大时的K值作为VMD的分解参数。

2）在确定K值后，选取峭度值大的模态分量作为将参与信号重构的模态。

3）将α从100~10000循环迭代寻找平均峭度达到最大时的α值作为VMD的分解参数。使用确定好的参数对原信号重新进行模态分解并完成信号重构。

4）使用凹凸型阈值小波变换对步骤3）中重构后的信号进行进一步的降噪，滤除残留的噪声。即完成了整个降噪流程。

2 PD信号降噪的仿真分析

2.1 PD信号的仿真

PD信号有单指数衰减、单指数振荡衰减、双指数衰减、双指数振荡衰减[17]4种数学模型，模型的表达式分别如下所示

式中：A1、A2为信号的2种幅值，本文分别设为2 MV、10 M V； τ为衰减系数，取5 μs；采样时间为50 μs；振荡频率fcl和fc2分别取5 MHz和10 MHz。在采样时间内仿真以上PD 4种信号，仿真信号及其频谱如图2所示。

图2 无噪声的PD仿真信号波形及频谱图Fig. 2 Simulate waveform and frequency spectrum of noisefree partial discharge signals

从图2中可以看到纯净PD仿真信号和其频率分布。红色部分主要是2个指数衰减信号的频率分布；而蓝色的部分则是2个指数振荡衰减信号的频率分布；fcl和fc2分别为这2个信号的频率峰值。

在现场测试时常见的噪声类型有周期窄带干扰和白噪声干扰，所以这里再仿真PD中添加（0，0.52）分布的高斯白噪声信号，以及周期窄带干扰信号，其表达式如下[18]

窄带的幅值Ai都选为0.1 MV，噪声频率fi分别设为1.2、8、16、20.7 MHz。添加2种噪声后的PD仿真信号如图3所示。

图3 含噪声的PD仿真信号波形及其频谱Fig. 3 Partial discharge simulate signal containing noise and its spectrum

从图3的波形中可以看到PD信号脉冲已经被噪声淹没；从频谱图中也可以看到白噪声信号分布在整个频带上，周期窄带干扰则分布在设定的频率位置。PD信号的频谱分布被淹没在其中，2个信号峰值已经无法分辨。

2.2 VMD分解及重构

针对染噪的PD仿真信号，VMD的K值、重构模态选择以及 α值的调试过程如下。

1）K值的调试。

对于仿真信号，本文在这里分别使用比较分解后的模态中心频率和求模态的最大平均峭度值2种调试方法来确定K值。这里对仿真的PD信号的VMD分解，K值分别取2、3、4、5、6、7，不同K值下各IMF函数的中心频率如表1所示，各模态的峭度值和平均峭度值如表2所示。

表2 不同K值下的IMF函数的峭度值Table 2 The kurtosis value of of intrinsic mode function under different values of K

在本文仿真的纯净PD信号中，2个指数衰减信号中心频率接近0.25 MHz，2个指数振荡衰减频率分别为5 MHz、10 MHz。从表1可以看到K=5时，IMF1、IMF2、IMF3 3个模态的中心频率和PD信号的频率非常接近。所以这里的K值选为5最为合适。

表1 不同K值下的IMF函数的中心频率Table 1 Center frequencies of intrinsic mode function under different values of K

从表2中可以看到,当K值为5时模态的平均峭度值最高。使用取模态平均峭度最大值时K值，K值同样应该为5。

2）参与重构的模态选择。

完成K值的调试后，这里根据峭度原则来进行对参与重构的模态筛选。K=5时各模态的峭度如表3所示。

从表3可以看出IMF1、IMF2、IMF3的峭度值远远大于剩余2个模态的峭度。这表明前3个模态中存在PD信号。因此选取IMF1、IMF2、IMF3作为将进行重构的模态。

表3 IMF的峭度值Table 3 The kurtosis value of intrinsic modes when K=5

3）惩罚因子 α值的调试。

在完成K值的选择和将参与重构的模态选择后，开始寻找合适的参数α值。本文使用一组仿真的PD信号进行参数α的调试，将α值从100~10000分别进行测试。测试中信噪比和平均峭度值变化如图4所示。

从图4中可以看到当惩罚因子α值达到5200时平均峭度值达到最大值，之后平均峭度值停止了增长。从6500开始出现下降，信噪比虽然一直是上升的状态，但α值在5200之后PD信号特征的损失一直是在增加的。因此这里α值取5200最为合适。

图4 惩罚因子与信噪比及峭度值的对应关系Fig. 4 Correspondence of penalty factor with signal-tonoise ratio and kurtosis value

在调整好参数后的IMF及IMF的频谱如图5所示，信号重构后的波形和频谱如图6所示。

从图5可以看到IMF1主要包含的是2个指数衰减信号。IMF2、IMF3 2个模态分量分别包含了2个指数振荡衰减信号；从图6可以看到，经过调整好参数的VMD分解重构后，大部分噪声已经被滤除出去了。频谱图和纯净的PD信号比较接近了，fcl和fc2两个频率峰值隐约可以辨认。但还残留一些和PD信号频带相同或相近的噪声信号。

图5 VMD分解IMF分量及频谱Fig. 5 VMD decomposition of IMF components and spectrogram

图6 重构后的信号及频谱Fig. 6 Reconstructed signals and spectrogram

2.3 凹凸型阈值的小波降噪

这里再通过凹凸型阈值小波对VMD分解重构后的PD信号进行进一步的降噪，进一步降噪后的波形和频谱如图7所示。

图7 进一步降噪后的波形及频谱图Fig. 7 Further noise reduction results and spectrogram

从图7可以看到经过凹凸型阈值小波的降噪处理后和PD信号频带相同或相近的残留噪声被进一步滤除。信号频率分布和纯净的局放信号更加接近，fcl和fc22个频率峰值更加清晰。

2.4 不同降噪方法降噪效果的对比

这里使用dB6小波硬阈值（方法1）、软阈值降噪（方法2）、VMD+硬阈值小波降噪（方法3）以及本文的降噪方法进行对比。本文选用的降噪评价指标分别为用信噪比(signal to noise ratio,SNR)、均方根误差(root-mean-square Error,RMSE)以及波形相似系数(normalized correlation coefficient, NCC)[19-20]。SNR参数值越高表明降噪效果越好；RMSE越接近0表明降噪后的信号畸变量越少；NCC越接近于1表明降噪后的波形和原始信号的波形越相似。各个降噪方法降噪效果的参数对比如表4所示。