自适应奇异值分解局放信号降噪方法

孙传铭， 魏 隆， 张梦楠， 刘 凯， 潘贵翔， 高国强

(1. 中车青岛四方机车车辆股份有限公司， 山东 青岛 266111；2. 西南交通大学电气工程学院， 四川 成都 610031)

1 引言

高压电气设备局部放电(Partial Discharge，PD)检测过程中往往会受到各种背景噪声的干扰，尤其在工程现场的电磁干扰对检测结果的影响更加严重，有时甚至会出现局放信号完全被背景噪声湮没的情况，对后续电气设备绝缘状态判断及检修带来一定的困难[1]。局放测试中的干扰主要分为以下三种[2-5]：随机性脉冲干扰、周期性窄带干扰和白噪声干扰。其中，随机脉冲干扰通常强度大、频率低，易于识别和滤除。周期性窄带干扰主要来源于电网络内部及环境中的无线电广播等信号，其出现频率高，幅值大，常在时域中将局放信号湮没，且在频域范围内经常与局放信号发生混叠，对局放信号的检测影响很大。白噪声干扰主要是电气设备的热噪声引起的宽带干扰随机信号，在频域上与局放信号具有相似特征。因此，如何有效滤除周期性窄带干扰和白噪声成为局放信号研究的一大难点。

针对周期性窄带干扰和白噪声干扰混合噪声的抑制方法，国内外学者进行了大量的研究。文献[6]提出结合广义S时频变换和奇异值分解去噪方法抑制信号中的混合噪声，该方法可有效地抑制混合噪声，但去除窄带干扰信号时需要人为判断窄带干扰区域，受人为因素影响存在一定误差，不具备自适应性。文献[7]针对奇异值分解耗时较长且有效奇异值数量难以选择问题，提出基于滑动短时数据能量窗的奇异值分解降噪方法，该方法无需预先假设信号中含有周期性窄带干扰，可自动实现周期性窄带干扰的甄别和混合噪声的抑制。但该方法在去除窄带干扰时往往受计算精度的影响而存在误差，最终影响降噪结果。文献[8]提出基于经验模态分解(Empirical Mode Decomposition，EMD)和独立成分分析(Independent Component Analysis，ICA)的高压电缆局放信号降噪方法，利用快速ICA方法进一步滤除含噪IMF分量中的噪声信号，但该方法并未明确提及是否具备同时滤除两种噪声的能力。文献[9]提出一种基于总体经验模态分解和补充总体经验模态分解的局部放电阈值去噪新方法，该方法通过对总体经验模态分解的IMF分量进一步通过补充总体经验模态分解提高降噪能力，但该方法耗时较长，应用受限。文献[10]提出基于经验小波和小波变换的局放信号降噪方法，通过两种方法结合实现降噪优化，但对于小波变换依旧存在基函数和分解层数的选择问题，自适应性能较差。

针对上述局放信号混合噪声干扰抑制存在的问题，本文提出一种自适应奇异值分解降噪方法。该方法首先对测试信号构建Hankel矩阵，以此作为轨迹矩阵进行奇异值分解[11]。通过提取前两个奇异值进行重构并结合功率谱熵自适应判断染噪信号中是否存在窄带干扰；随后确定奇异值系列特征量，结合K类均值聚类[12]对窄带干扰、有效PD信号和白噪声所对应的奇异值进行划分，对有效PD信号对应的奇异值进行重构进而还原PD信号。该方法可自适应地判断是否存在窄带干扰，从而决定聚类区间；通过1次K类均值聚类分类即可获取有效奇异值数据，自适应性能良好。

2 奇异值分解

2.1 奇异值分解原理

2.1.1 轨迹矩阵的构建

本文选取Hankel矩阵作为奇异值分解的轨迹矩阵。Hankel矩阵具体构建方式如下：

设染噪信号X为：

X=[x(1),x(2),…，x(N)]

(1)

对采样系列X构造Hankel矩阵：

(2)

式中，N=m+n-1，本文中n取为N/2。

2.1.2 奇异值获取

矩阵A是一个m×n的矩阵，其秩为r，则必存在m×m的正交矩阵U和n×n的正交矩阵V，使得：

A=UΛVT

(3)

其中

(4)

对角矩阵Λ除了前r阶对角元素外，其他元素均为零。对角元素∂i即为矩阵A的奇异值，且数值由大到小排列，奇异值的大小反映了能量的集中情况。通过对窄带干扰、有效PD信号和白噪声对应的奇异值规律进行剖析，进而选取合适的奇异值进行重构，即可还原真实的局放信号。

2.2 信号仿真

结合以往数据仿真经验，本文选取单指数振荡衰减模型和双指数振荡衰减模型来模拟理想局放信号[13]。具体表达式如下：

(5)

(6)

式中，f为振荡频率；τ为衰减系数；A1、A2为脉冲幅值。

本文仿真了四种局部放电脉冲，其中脉冲模型1和模型2根据式(5)得出，脉冲模型3和模型4根据式(6)得出。脉冲仿真参数见表1。

表1 局放仿真信号参数Tab.1 PD simulation signal parameters

周期性窄带干扰通常呈正弦或余弦波形[14]，且PD信号实际检测中往往会存在窄带干扰与PD信号混叠的问题。故本文选取窄带干扰频率分别0.5 MHz、1 MHz、2 MHz、5 MHz、7 MHz。周期性窄带干扰的具体表达式如下：

(7)

式中，Ai对应各窄带干扰信号幅值；fi为频率。模拟窄带干扰仿真信号参数见表2。

表2 窄带干扰仿真信号参数Tab.2 Parameters of narrow-band interference simulation signals

实际运行环境中除受窄带干扰影响外，往往还会受到白噪声的干扰，白噪声利用高斯白噪声模拟产生。仿真获取理想PD仿真信号如图1(a)所示，添加周期性窄带干扰和白噪声后的信号如图1(b)所示，图1(c)为染噪信号频域谱图。

图1 仿真波形Fig.1 Simulation waveforms

2.3 窄带干扰奇异值特征分析

针对窄带干扰与奇异值分解间存在的关系问题，文献[15]发现对含有单个频率窄带干扰的PD信号进行SVD分解时，提取前两个奇异值可有效提取窄带干扰信号。文献[16]发现每个频率的窄带干扰都对应两个非0奇异值，通过提取窄带干扰频率个数n对应的前2n个奇异值即可提取窄带干扰并滤除。文献[17]进一步指出随着采样数据长度增加，窄带干扰对应奇异值幅值越来越大，而局放信号对应奇异值变化较小。因此，通过增加数据长度，可保证窄带干扰被全部滤。染噪信号奇异值随数据长度变化情况如图2所示。

图2 奇异值与数据长度的关系Fig.2 Relation between singular values and data length

3 局放混合噪声抑制方法

3.1 窄带干扰判别

为实现对局放信号混合噪声的自适应抑制，首先需要对染噪信号中是否存在窄带干扰进行判别。根据前人对窄带干扰与PD信号奇异值规律的剖析，本文提取信号前两个奇异值进行重构，根据重构后的信号是否符合正(余)弦规律即可判断是否存在窄带干扰。对于功率谱熵而言，信号混乱程度越高，其功率谱熵越大，混乱程度越低，功率谱熵越小。正弦信号混乱程度较局放信号较小，因此，本文引入功率谱熵[18,19]的概念对正(余)弦信号进行检测。利用正弦信号与局放信号间功率谱熵大小的差异，判断是否存在窄带干扰。基于功率谱熵检测的具体步骤为：

(1)将信号x(t)经FFT变换得到功率谱为：

(8)

(9)

式中，N为数据点个数。

(2)求取信号功率谱熵H为：

(10)

(11)

(3)H作为检测统计量为：

(12)

当检测统计量H小于检测阈值T时，即可确定重构信号是窄带干扰。

图3为局放信号和窄带干扰信号幅值和频率改变时分别对应的功率谱熵幅值。从图3中可以看出，无论窄带干扰信号幅值和频率如何变化，对应功率谱熵幅值均小于1；而局放信号功率谱熵幅值始终大于1。经过多组数据分析，本文最终设定检测阈值T=1，若重构信号对应检测阈值T<1，判定染噪信号中含有窄带干扰。

图3 功率谱熵求取结果Fig.3 Result of power spectrum entropy obtained

3.2 有效奇异值选取

奇异值有效个数选取问题一直是奇异值分解降噪的关键[20]。如果有效奇异值个数选取过少，将会损失局放信号部分有用信息；如果有效奇异值个数选取过多，则降噪效果不明显。此外，如何实现有效奇异值个数的自适应选取也是一大研究重点。根据以往的研究得知，窄带干扰信号奇异值远大于局放混合白噪声信号对应奇异值，且该数值位于奇异值分解前列。同时，局放信号相对白噪声而言，其奇异值相对较大且数据较分散。因此，本文引入K类均值聚类算法[21]对局放混合噪声进行抑制。为凸显奇异值大小，本文选取奇异值本身F1作为奇异值系列特征参量，同时，为表征数据离散程度，同时选取奇异值子集标准差F2作为另一奇异值系列特征量，并与文献[17]中选取以奇异熵增量F3及其能量F4作为奇异值系列特征量进行比较，对混合染噪信号进行分类，结果如图4和图5所示。各类特征参量计算公式及表达式如下。

图4 F1、F2特征参量K类均值聚类结果Fig.4 Results of K-means clustering of characteristic parameters F1 and F2

图5 F3、F4特征参量K类均值聚类结果Fig.5 Results of K-means clustering of characteristic parameters F3 and F4

奇异值子集ssi构建及子集标准偏差σi计算公式为：

ssi=[∂r,…,∂i]

(13)

(14)

式中，∂i对应第i个奇异值；μi为ssi数据均值。

奇异熵增量计算公式：

(15)

式中，r为奇异值总数。

奇异熵增量能量计算公式：

ei2=(Δei)2

(16)

特征量F1、F2、F3、F4表达式如下：

F1=[∂r,∂r-1,…,∂1]

(17)

F2=[σ1,σ2,…,σr]

(18)

F3=[Δe1,Δe2,…,Δer]

(19)

F4=[Δe12,Δe22,…,Δer2]

(20)

从图4中可以看出，本文选取的奇异值特征量实现了窄带干扰、有效PD信号和白噪声信号对应奇异值的有效分类，重构信号在保证滤除混合噪声的同时保留了局放信号的完整性，本文方法更适用于同时实现3种不同信号的有效分离。此外，在局放信号中混叠窄带干扰的情况下，本文分别比较了2次K类均值聚类分类次数为2的奇异值分类和1次K类均值聚类分类次数为3的分类。结果表明，采用本文方法选取的特征参量进行两种分类算法获取的有效奇异值基本吻合。为节省计算时间，本文最终选取仅作1次K类均值聚类分类次数为3的计算。

3.3 混合噪声抑制步骤

要实现混合噪声的自适应抑制，窄带干扰信号的判别至关重要。本文首先对信号进行奇异值分解，通过提取前两个奇异值重构判断信号中是否存在窄带干扰；然后采用K类均值聚类最终提取有效奇异值，进而获取局放信号。混合噪声抑制步骤具体如下：

(1)对染噪信号构建Hankel矩阵，以此作为轨迹矩阵进行奇异值分解；

(2)重构前两个奇异值获取重构信号，判断信号对应功率谱熵T是否大于1。若T<1，则存在窄带干扰，设置K类均值聚类分类个数n=3；若T>1，则不存在窄带干扰，设置n=2；

(3)分别以奇异值本身F1和奇异值子集标准差F2为奇异值系列特征量，通过K类均值聚类将奇异值系列分为n类；

(4)若n=3，则选取第二类奇异值数据进行重构还原PD信号；若n=2，则选取第一类奇异值数据进行重构还原PD信号。

综上所述，本文还原PD信号的流程图如图6所示。

图6 本文降噪方法流程图Fig.6 Flow chart of noise reduction method in this paper

4 去噪效果对比

为分析本文的自适应奇异值分解降噪方法对局放信号的降噪效果，对原始PD仿真信号加入周期性窄带干扰并叠加分布为(0，10)的高斯白噪声进行降噪处理。通过引入FFT-小波变换降噪、S时频变换-EEMD联合去噪方法与本文方法进行对比。各方法降噪结果如图7所示。经对比得出，FFT-小波变换降噪可有效滤除噪声，但局放信号也被部分滤除，导致信号减小；S变换-EEMD联合降噪既不能保证噪声的高精度滤除，同时PD信号存在部分衰减。本文所选方法能同时满足噪声信号的高精度滤除和PD信号的高度还原。

图7 3种方法降噪结果对比Fig.7 Denoising results of three method

本文引入去噪评价参数信噪比(Signal to Noise Ratio，SNR)、均方误差(Mean Square Error，MSE)和波形相似参数(Normalized Correlation Coefficient，NCC)[22]进一步对降噪效果进行评估。

去噪评价参数计算结果如表3所示。从表3中可以看出，本文所采用的方法具有明显的优势，无论从信噪比、均方误差还是波形相似参数上都显示出非常好的效果，对于混合噪声干扰的抑制效果最好，且信号还原度最高。

表3 去噪评价参数计算结果Tab.3 Calculation results of denoising evaluation parameters

5 实测信号去噪分析

为检验本文方法对于实测PD信号滤除混合噪声的能力。基于实验室条件下搭建电缆终端刀痕缺陷测试模型如图8所示，测试采用的高频脉冲电流传感器-6 dB带宽为80 kHz～40 MHz，采样率为50 MSa/s。测试得局放波形如图9(a)所示。实测PD信号基于理想试验条件下测得，故而PD信号明显，而环境噪声干扰很小。为测试本文降噪方法对实测信号的去噪效果，通过对实测PD信号施加3个幅值为0.5 mV，频率分别为0.5 MHz、2 MHz和8 MHz的周期性窄带干扰信号，并叠加分布为(0，10)的高斯白噪声。染噪信号如图9(b)所示。

图8 局放检测平台原理Fig.8 Schematic of PD detection circuit

图9 实验室实测PD信号Fig.9 Measured PD signals in laboratory

分别采用FFT-小波变换降噪、S时频变换-EEMD联合去噪方法和本文方法对添加周期性窄带干扰和白噪声的实验室环境下实测PD信号进行降噪处理，各方法降噪结果如图10所示。从图10中可以很明显地看出，本文降噪方法能够高度还原实测PD信号，且抑制噪声效果最好。FFT-小波变换降噪明显改变了PD信号的特征。同时，PD信号明显减小；S时频变换-EEMD联合降噪方法较FFT-小波变换降噪效果更佳，但同样存在信号衰减问题。此外，S时频变换-EEMD降噪去除噪声效果相对较差。

图10 实验室实测PD信号降噪结果Fig.10 Noise reduction results of PD signal measured in laboratory

由于无法测得完全不含噪声的PD信号，故无法使用上述去噪评价参数对各方法降噪效果进行定量分析。因此，本文引入噪声抑制比μ1和幅值衰减比μ2对降噪效果进行评价[11]。其中，μ1反映了降噪后信号的凸显程度。μ1越大，说明降噪方法去噪效果越好。μ2反映了降噪前后PD信号的衰减程度。μ2越大，说明降噪后PD信号衰减越严重。μ1、μ2具体定义如式(21)、式(22)所示。

μ1=10(lgδ12-lgδ22)

(21)

(22)

式中，δ1、δ2分别为降噪前后信号的标准偏差;Am1、Am2分别为降噪前后信号的最大幅值。各降噪方法降噪评价参数计算结果如表4所示。

从表4中可以看出，文本降噪方法降噪效果最好，信号衰减程度最小。故而选取本文的方法在降噪上占有很大的优势。

表4 降噪评价参数计算结果Fig.4 Results of evaluation parameters of noise reduction

6 结论

本文基于奇异值分解自适应降噪，提出了一种有效滤除白噪声和周期性窄带干扰的降噪方法。通过与FFT-小波变换降噪和S时频变换-EEMD联合降噪方法进行对比，分析结果发现本文方法具有更优的降噪效果，且PD信号衰减最小，还原度最高。具体结论如下：

(1)信号降噪前首先对染噪信号提取前两个奇异值重构，自适应判断信号中是否含有窄带干扰，避免了因窄带干扰存在与否问题导致的降噪失误，从而为后续利用SVD实现降噪奠定基础。

(2)利用Hankel矩阵作为轨迹矩阵进行奇异值分解时，窄带干扰频率个数对应两倍奇异值个数，且采样数据足够长时，窄带干扰对应奇异值数值远大于PD信号。

(3)采用奇异值本身和奇异值子集标准差作为奇异值系列特征量时，利用K类均值聚类分类方法可一次实现对窄带干扰、有效PD信号和白噪声对应奇异值的有效分离，且分类效果与分别进行两次分类的效果一致，分类所需用时有效缩短。

(4)根据窄带干扰判断结果可自适应确定K类均值分类个数并进行分类，进而对有效奇异值进行重构得到PD信号，自适应性能良好。

(5)通过与FFT-小波变换降噪和S时频变换-EEMD联合降噪方法进行对比，发现本文方法能更好地抑制噪声，同时保证PD信号的高度还原。