CN113887362A - 一种局部放电信号的特征提取方法 - Google Patents

Abstract

本发明涉及一种局部放电信号的特征提取方法，包括以下步骤：局部放电信号进行采集；根据采集到的局部放电信号构建Hankel矩阵，并进行奇异值分解，获取重构后的轨迹矩阵；从原始局部放电信号中减去周期性窄带干扰，获取初步提取局部放电信号；将初步提取局部放电信号分解为多个IMF分量，利用相关系数滤除白噪声；计算保留的IMF分量的样本熵，将样本熵的中值作为阈值，拾取出保留的IMF分量中的高频分量，对拾取出的高频分量进行去噪，将经过改进阈值去噪的高频分量与余下的低频分量重构，得到最终提取的PD信号。与现有技术相比，本发明结合硬阈值函数和软阈值函数对IMF分量进行阈值处理，提高特征提取效果。

Description

一种局部放电信号的特征提取方法

技术领域

本发明涉及电气工程在线监测领域，尤其是涉及一种局部放电信号的特征提取方法。

背景技术

局部放电(Partial Discharge,简称PD)是导致变压器等电力设备绝缘老化或损坏的主要原因，对电力设备进行局部放电的在线监测是保证电网安全、稳定、连续运行的关键手段。但现场的局部放电信号中往往包含着强烈的噪声干扰，其中又以周期性窄带干扰和高斯白噪声最为严重。局部放电初期产生的信号较弱，导致有效的局部放电信号淹没在强烈的背景干扰中，致使局部放电信号难以识别。因此，对PD信号进行干扰抑制的研究具有重要意义。

现有技术中常用数字信号处理的方法抑制PD信号的干扰，主要有：广义形态学滤波法、小波分解法、经验模态分解法(EMD)和逆向分离法。广义形态学滤波法虽然可以有效抑制窄带干扰，但存在结构元素的选取问题；小波分解法通过时频域描述信号的特征，可以显著抑制白噪声，但是其性能依赖于小波函数、去噪阈值和分解层数的选取，缺乏自适应性；经验模态分解法(EMD)可以自适应的将信号分解为不同频率的本征模态函数(IMF)，根据局部放电信号与噪声信号在IMF上分布不同实现抑制干扰的效果，但是容易出现模态混叠的现象，在低信噪比下性能不稳定；逆向分离法结合独立分量分析进行盲源分离，但是其幅值存在较大的畸变。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种局部放电信号的特征提取方法。

本发明的目的可以通过以下技术方案来实现：

一种局部放电信号的特征提取方法，包括以下步骤：

S1：获取局部放电信号；

S2：根据采集到的局部放电信号构建Hankel矩阵，并进行奇异值分解，获取重构后的轨迹矩阵；

S3：根据重构后的轨迹矩阵分离出周期性窄带干扰，从原始局部放电信号中减去周期性窄带干扰，获取初步提取局部放电信号；

S4：利用CEEMDAN算法将初步提取局部放电信号分解为多个IMF分量，计算IMF分量与分解前信号的相关系数，当相关系数第一次出现局部极大值时剔除小于该阶的IMF分量；

S5：计算保留的IMF分量的样本熵，将样本熵的中值作为阈值，拾取出保留的IMF分量中的高频分量，对拾取出的高频分量进行基于Sigmoid函数的改进阈值去噪，将经过改进阈值去噪的高频分量与余下的低频分量重构，得到最终提取的局部放电信号。

优选地，所述步骤S1中采集的局部放电信号为：

X_i＝{x_i(1)，x_i(2)，…，x_i(N)}

其中，x_i(N)为局部放电信号的采样序列中第N个采样点，N为采样序列中采样点数量，

步骤S2中的Hankel矩阵为：

所述Hankel矩阵为m行k列矩阵，其中，m＝N+1-w，w＝N/2。

优选地，所述重构后的轨迹矩阵为H_s：

其中，p为表征周期性窄带干扰的奇异值点个数，U_p为左奇异矩阵的前p列，S_p为前p个奇异值构成的对角矩阵，

为右奇异矩阵的前p行。

优选地，所述的步骤S3具体包括：

S31：将重构后的轨迹矩阵H_s展开：

其中，y_s为轨迹矩阵H_s的元素

S32：求取反对角元素的均值对矩阵降维得到一维时间序列：

Y_i＝{y_i(1)，y_i(2)，…y_i(N)}

其中，

S33：将重构得到的一维时间序列从原始局部放电信号中减去得到抑制了窄带干扰的初步提取局部放电信号Z_i：

Z_i＝X_i-Y_i。

其中，X_i为局部放电信号。

优选地，所述利用CEEMDAN算法将初步提取局部放电信号分解为多个IMF分量的具体步骤包括：

在初步提取局部放电信号Z_i中加入白噪声ωⁱ(t)，得到新的信号Z_i(t)+ε₀ωⁱ(t)，ε₀为噪声系数，对新的信号用EMD方法进行I次分解，求取均值得到第一个IMF分量：

一阶残余分量为r₁(t)：r₁(t)＝Z_i(t)-IMF₁(t)

对信号r₁(t)+ε₁E₁(ωⁱ(t))再重复进行I次分解，其中E_j(·)为信号经过EMD分解后的第j个固有模态函数，得到第二个IMF分量：

对于k＝2，…，K，计算k阶残余量为：r_k(t)＝r_k-1(t)-IMF_k(t)，提取信号r_k(t)+ε_kE_k(ωⁱ(t))中第一个IMF的分量，第k+1个模态函数分量定义为：

依次获取各阶IMF分量直到信号为单调函数，得到K个模态函数时，残差表示为：

初步提取局部放电信号Z_i表示为：

优选地，所述的相关系数的计算公式为：

其中，ρ为相关系数，

为初步提取局部放电信号Z_i的均值，

为第k阶IMF分量的均值。

优选地，所述步骤S5中将经过改进阈值去噪的高频分量与余下的低频分量重构的具体公式为：

其中，IMF′_k(t)为重构分量，T_k为阈值，S(α)为Sigmoid函数，

优选地，所述的阈值T_k为：

其中，C为常数；N为IMF序列的长度，E_k为样本熵提取的高频分量中第k阶IMF的能量。

优选地，所述的E_k的计算公式为：

其中，E^*为样本熵提取的高频分量中阶数最低的IMF对应的能量。

优选地，所述的E^*的计算公式为：

其中，IMF^*对应样本熵提取的高频分量中阶数最低的本征模态函数。

与现有技术相比，本发明具有如下优点：

(1)本发明结合Hankel矩阵构建、奇异值分解，并对信号进行自适应噪声的完备经验模态分解，能够有效对局部放电信号的噪声进行去除，自适应性强，在低信噪比状态下性能稳定，提取到的特征能更好的表征局部放电信号；

(2)本发明利用Hankel矩阵构建、奇异值分解，通过奇异值的局部极大值点确认周期性窄带干扰和PD信号的奇异值分界点，进而去除局部放电信号中的周期性窄带干扰，对局部放电信号进行初步有效去噪，为后续基于自适应噪声的完备经验模态分解提供更适用的局部放电信号。

(3)本发明利用CEEMDAN算法将初步提取局部放电信号分解为多个IMF分量，并计算相关系数，通过相关系数确定噪声和局部放电信号的分解，滤除大部分白噪声，最后计算样本熵对高频分量进行改进阈值滤波去除残余白噪声，最后与低频分量重构获取最终局部放电信号，对初步提取局部放电信号中的白噪声进行有效过滤，提高特征提取效果。

附图说明

图1为本发明局部放电信号特征提取的流程图。

图2为本发明实施例中的仿真局部放电信号示意图。

图3为本发明实施例中染噪局部放电信号示意图。

图4为对图3信号进行奇异值分解后奇异值的变化曲线示意图。

图5为图3信号经过Hankel-SVD处理后的信号示意图。

图6为图5信号经过CEEMDAN分解后形成的各个IMF分量的三维示意图。

图7为图6中各个IMF分量的时域波形及对应频谱示意图。

图8为最终提取的仿真局部放电信号示意图。

图9为实测局部放电信号示意图。

图10为对图9信号添加噪声后的染噪信号示意图。

图11为实际局部放电信号和去噪后局部放电信号对比的示意图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。注意，以下的实施方式的说明只是实质上的例示，本发明并不意在对其适用物或其用途进行限定，且本发明并不限定于以下的实施方式。

实施例

一种局部放电信号的特征提取方法，如图1所示，包括以下步骤：

S1：局部放电信号进行采集；

S2：根据采集到的局部放电信号构建Hankel矩阵，并进行奇异值分解，获取重构后的轨迹矩阵，具体为：

构建Hankel矩阵，采样序列为：X_i＝{x_i(1)，x_i(2)，…x_i(N)}，其中N为序列长度，i为input的缩写，表示为输入信号。将PD信号的采样序列构建为Hankel矩阵，构建方式为：

其中：m＝N+1-w，w＝N/2。

奇异值分解的原理为：

设有双线性函数：f(x，y)＝x^TAy

其中，x，y∈R^n×1，A∈R^n×1，引入线性变化：x＝Uξ，y＝Vη，U，V∈R^n×n，上式变化为：

f(x，y)＝(Uξ)^TA(Vη)＝ξ^TU^TAVη

令s＝U^TAV，有：

f(x，y)＝ξ^TSη

通过选择合适的正交矩阵U、V，可以使得S为正交矩阵，即：

S＝diag(σ₁，σ₂，…，σ_n)

公式两边左乘U、右乘V^T得到：

A＝USV^T

上式被称为矩阵A的奇异值分解。奇异值分解(SVD)是线性代数中的一个重要概念，SVD可以将任意维实矩阵分解为左奇异矩阵U、对角矩阵S和右奇异矩阵V^T的乘积，其中对角矩阵S主对角线上的元素称为奇异值。

周期性窄带干扰不仅能量远高于PD信号，而且其频谱与PD信号存在重叠部分，难以得到抑制。奇异值分解(SVD)利用周期性窄带干扰与PD信号正交不相关的特点，将原始PD信号分解为不同的正交分量。表征窄带干扰的奇异值因为携带大量能量会显著大于其他奇异值。奇异值分解的对象是矩阵，不可直接对原始PD信号采样序列进行分解。利用Hankel矩阵将一维PD信号采样序列构建至二维的矩阵空间，再对Hankel矩阵进行SVD分解即可获得PD信号对应的左右奇异矩阵和奇异值。

为了找到有效奇异值信号进行重构，计算奇异值变化的斜率g_n：

其中：(p_n，q_n)是第n个奇异值的坐标，(p_n+1，q_n+1)是第n+1个奇异值的坐标。将斜率最大值点作为分界点，保留表征周期性窄带干扰的奇异值点而将其余奇异值置零进行重构，设重构后得到的轨迹矩阵为H_s：

其中，p为表征周期性窄带干扰的奇异值点个数，U_p为左奇异矩阵的前p列，S_p为前p个奇异值构成的对角矩阵，

为右奇异矩阵的前p行。

S3：根据重构后的轨迹矩阵分离出周期性窄带干扰，从原始局部放电信号中减去周期性窄带干扰，获取初步提取局部放电信号：

将轨迹矩阵展开：

其中，y_s为轨迹矩阵的元素。求取反对角元素的均值对矩阵降维得到一维时间序列Y_i＝{y_i(1)，y_i(2)，…y_i(N)}。其中：

将重构得到的时间序列Y_i从原始采样序列X_i中减去得到抑制了窄带干扰的初步提取局部放电信号信号Z_i：

Z_i＝X_i-Y_i

S4：利用CEEMDAN算法将初步提取局部放电信号分解为多个IMF分量，计算IMF分量与分解前信号的相关系数，当相关系数第一次出现局部极大值时剔除小于该阶的IMF分量，具体为：

CEEMDAN的原理为：

在初步提取局部放电信号Z_i中加入不同的白噪声ωⁱ(t)，得到新的信号Z_i(t)+ε₀ωⁱ(t)，ε₀为噪声系数。对新的信号用EMD方法重复进行I次分解，求取均值得到第一个固有模态函数(Instrinsic mode function，IMF)：

第一个的残余分量则为：

r₁(t)＝Z_i(t)-IMF₁(t)

定义算子E_j(·)为信号经过EMD分解后的第j个固有模态函数。对信号r₁(t)+ε₁E₁(ωⁱ(t))再重复进行I次分解，得到第二个IMF分量：

对于k＝2，…，K，计算k阶残余量为：r_k(t)＝r_k-1(t)-IMF_k(t)

提取信号r_k(t)+ε_kE_k(ωⁱ(t))中第一个IMF的分量，第k+1个模态函数分量定义为：

以此类推得到各阶IMF分量，直到信号为单调函数时停止分解，当得到K个模态函数时，残差表示为：

原始信号x(t)表示为：

本发明选择当相关系数第一次出现局部极大值时，将该阶IMF分量作为噪声分量和PD信号分量的分界，之后将噪声分量剔除实现对白噪声的抑制。相关系数的计算公式为：

其中，ρ为相关系数，

为初步提取局部放电信号Z_i的均值，

为第k阶IMF分量的均值。

S5：计算保留的IMF分量的样本熵，将样本熵的中值作为阈值，拾取出保留的IMF分量中的高频分量，对拾取出的高频分量进行基于Sigmoid函数的改进阈值去噪，将经过改进阈值去噪的高频分量与余下的低频分量重构，得到最终提取的PD信号。

本发明计算经过相关系数法舍弃掉噪声分量后的各个IMF分量的样本熵，求取样本熵的中值作为划分高频分量和低频分量的阈值。样本熵大于阈值的是序列复杂性高的高频分量，样本熵低于阈值的是序列复杂性低的低频分量。样本熵的计算步骤如下：

假设长度M的序列x(t)＝{x₁，x₂，…x_M}，参照原始序列重构得到向量X(i)＝{x_i，x_i+1，…，x_i+u-1}，其中1≤i≤M-u，u为嵌入维数；

计算d_ij(1≤j≤M-u，j≠i)，d_ij为X(i)和X(j)对应元素差值绝对值的最大值；

统计d_ij中小于r的个数记为num(d_ij＜r)，r为相似容限。定义B_iu(r)＝num(d_ij＜r)/(M-u-1)；

求B_iu(r)的平均值记为B^u(r)；

对维数u+1重复上述步骤得到B_iu+1(r)和B^u+1(r)；

计算样本熵：

受小波阈值法启发，对CEEMDAN分解得到的IMF分量可以做类似阈值处理，有两种形式，硬阈值函数：

软阈值函数：

本发明针对PD信号的特征，结合硬阈值函数和软阈值函数的优缺点，提出一种改进的阈值函数用于CEEMDAN分解后对IMF分量的阈值处理：

其中：

对阈值T_k采用下述模型进行自适应选取：

其中C为常数；N为IMF序列的长度，E_k为样本熵提取的高频分量中第k阶IMF的能量，可以使用下式进行估算：

采用低阶IMF的能量估计剩余IMF的能量，E^*为样本熵提取的高频分量中阶数最低的IMF对应的能量：

其中IMF^*对应样本熵提取的高频分量中阶数最低的本征模态函数。

本实施例中，局部放电信号为振荡衰减信号，采用单指数衰减振荡模型和双指数衰减振荡模型进行描述。

单指数衰减振荡型：

f(t)＝Ae^-t/τsin(2πf_ct)

双指数衰减振荡型：

A(e^-1.3t/τ-e^-2.2t/τ)sin(2πf_ct)

其中，A为幅值；τ为衰减系数，f_c为振荡频率。采样频率为20MHz，PD脉冲1、3为单指数衰减振荡形式，PD脉冲2、4为双指数衰减振荡形式。图2中PD信号的仿真参数如表1所示：

表1 PD信号仿真参数

PD脉冲标记	1	2	3	4
					振荡频率f<sub>c</sub>/MHz	1	1.5	1	1.5
信号幅值A/mV	0.3	1.3	0.3	1.3
					衰减系数τ/us	2	3	2	3

周期性窄带干扰主要由载波通信引起，由不同频率的正弦信号叠加形成，数学表达式为：

其中，h为周期性窄带干扰的数量，A_i为周期性窄带干扰幅值，f_i为周期性窄带干扰的频率。仿真选取20MHz采样频率，在叠加白噪声基础上再叠加五种窄带干扰，频率为：100kHz、500kHz、1.2MHz、5MHz和7MHz。图3为最终染噪信号，此时PD脉冲已经完全淹没在周期性窄带干扰和白噪声中。

图4为对图3进行Hankel-SVD处理后的奇异值变化曲线，斜率突变点发生在第十个奇异值点处，即前十个奇异值队应的即是周期性窄带干扰，对表征周期性窄带干扰的奇异值进行重构得到轨迹矩阵。对重构得到的轨迹矩阵反对角线求均值分离出周期性窄带干扰的一维时间序列，从原始PD信号中减去重构的周期性窄带干扰得到初步提取的PD信号，初步提取结果如图4所示。

图6为图5信号经过CEEMDAN分解后形成的前十个IMF分量的三维示意图，可以看到PD脉冲和白噪声初步得到了分离。图7为图6中各个IMF分量的时域波形及对应频谱示意图，前三个IMF分量频谱范围较宽，对应着分离出的大部分白噪声。各个IMF分量的相关系数及样本熵计算结果如表2所示，相关系数在IMF3时首次出现局部极大值，剔除掉IMF1、IMF2分量；利用样本熵划分后，需要进行改进阈值去噪的高频分量为：IMF3、IMF4、IMF5、IMF6、IMF8，余下低频分量为：IMF7、IMF9、IMF10、IMF11、IMF12。

表2 相关系数及样本熵计算结果

固有模态分量	相关系数	样本熵
			IMF1	0.0024	2.1309
IMF2	0.0212	2.1705
			IMF3	0.4408	1.2262
IMF4	0.4341	1.3283
			IMF5	0.6465	1.1693
IMF6	0.6370	0.5545
			IMF7	0.3878	0.4188
IMF8	0.0480	0.5338
			IMF9	0.0070	0.1714
IMF10	-0.0224	0.0717
			IMF11	0.0120	0.0379
IMF12	-0.0038	0.0033

图8为最终提取的仿真局部放电信号，其中的周期性窄带干扰和白噪声已经得到了有效去除。

为了证明本发明的有效性，结合附图给出实际局部放电信号的处理结果。

图9为实测局部放电信号，图10为对图9信号添加噪声后的染噪信号。图11为实际局部放电信号和去噪后局部放电信号对比，可以看到本发明从含有高强度噪声的染噪信号中提取得到的PD信号与实测PD信号的波形高度重合，畸变较小，在波形细节和上升下降趋势上也有着很高的一致性。

因为无法获取完全不受噪声干扰的实测PD信号，引用噪声抑制比对去噪性能进行评估，噪声抑制比定义为：

其中，σ₁为去噪前信号的标准差，σ₂为去噪后信号的标准差。噪声抑制比反映了去噪前后PD信号的凸显程度，其值越大越好。各种方法对实测信号处理后的噪声抑制比如表3所示，对比可知，本发明的噪声抑制比明显优于其他方法。

表3 实测信号指标对比

指标	CEEMD-小波阈值	EEMD-SVD	本文方法
				噪声抑制比	17.6355	22.8654	29.4302

上述实施方式仅为例举，不表示对本发明范围的限定。这些实施方式还能以其它各种方式来实施，且能在不脱离本发明技术思想的范围内作各种省略、置换、变更。

Claims (10)

1.一种局部放电信号的特征提取方法，其特征在于，包括以下步骤：

S1：获取局部放电信号；

S2：根据采集到的局部放电信号构建Hankel矩阵，并进行奇异值分解，获取重构后的轨迹矩阵；

S3：根据重构后的轨迹矩阵分离出周期性窄带干扰，从原始局部放电信号中减去周期性窄带干扰，获取初步提取局部放电信号；

S4：利用CEEMDAN算法将初步提取局部放电信号分解为多个IMF分量，计算IMF分量与分解前信号的相关系数，当相关系数第一次出现局部极大值时剔除小于该阶的IMF分量；

S5：计算保留的IMF分量的样本熵，将样本熵的中值作为阈值，拾取出保留的IMF分量中的高频分量，对拾取出的高频分量进行基于Sigmoid函数的改进阈值去噪，将经过改进阈值去噪的高频分量与余下的低频分量重构，得到最终提取的局部放电信号。

2.根据权利要求1所述的一种局部放电信号的特征提取方法，其特征在于，所述步骤S1中采集的局部放电信号为：

X_i＝{x_i(1),x_i(2),…,x_i(N)}

其中，x_i(N)为局部放电信号的采样序列中第N个采样点，N为采样序列中采样点数量，

步骤S2中的Hankel矩阵为：

所述Hankel矩阵为m行k列矩阵，其中，m＝N+1-w，w＝N/2。

3.根据权利要求1所述的一种局部放电信号的特征提取方法，其特征在于，所述重构后的轨迹矩阵为H_s：

其中，p为表征周期性窄带干扰的奇异值点个数，U_p为左奇异矩阵的前p列，S_p为前p个奇异值构成的对角矩阵，

为右奇异矩阵的前p行。

4.根据权利要求1所述的一种局部放电信号的特征提取方法，其特征在于，所述的步骤S3具体包括：

S31：将重构后的轨迹矩阵H_s展开：

其中，y_s为轨迹矩阵H_s的元素

S32：求取反对角元素的均值对矩阵降维得到一维时间序列：

Y_i＝{y_i(1),y_i(2),…y_i(N)}

其中，

S33：将重构得到的一维时间序列从原始局部放电信号中减去得到抑制了窄带干扰的初步提取局部放电信号Z_i：

Z_i＝X_i-Y_i。

其中，X_i为局部放电信号。

5.根据权利要求1所述的一种局部放电信号的特征提取方法，其特征在于，所述利用CEEMDAN算法将初步提取局部放电信号分解为多个IMF分量的具体步骤包括：

在初步提取局部放电信号Z_i中加入白噪声ωⁱ(t)，得到新的信号Z_i(t)+ε₀ωⁱ(t)，ε₀为噪声系数，对新的信号用EMD方法进行I次分解，求取均值得到第一个IMF分量：

一阶残余分量为r₁(t)：r₁(t)＝Z_i(t)-IMF₁(t)

对信号r₁(t)+ε₁E₁(ωⁱ(t))再重复进行I次分解，其中E_j(·)为信号经过EMD分解后的第j个固有模态函数，得到第二个IMF分量：

对于k＝2,…,K，计算k阶残余量为：r_k(t)＝r_k-1(t)-IMF_k(t)，提取信号r_k(t)+ε_kE_k(ωⁱ(t))中第一个IMF的分量，第k+1个模态函数分量定义为：

依次获取各阶IMF分量直到信号为单调函数，得到K个模态函数时，残差表示为：

初步提取局部放电信号Z_i表示为：

6.根据权利要求5所述的一种局部放电信号的特征提取方法，其特征在于，所述的相关系数的计算公式为：

其中，ρ为相关系数，

为初步提取局部放电信号Z_i的均值，

为第k阶IMF分量的均值。

7.根据权利要求5所述的一种局部放电信号的特征提取方法，其特征在于，所述步骤S5中将经过改进阈值去噪的高频分量与余下的低频分量重构的具体公式为：

其中，IMF′_k(t)为重构分量，T_k为阈值，S(α)为Sigmoid函数，

8.根据权利要求7所述的一种局部放电信号的特征提取方法，其特征在于，所述的阈值T_k为：

其中，C为常数；N为IMF序列的长度，E_k为样本熵提取的高频分量中第k阶IMF的能量。

9.根据权利要求8所述的一种局部放电信号的特征提取方法，其特征在于，所述的E_k的计算公式为：

其中，E^*为样本熵提取的高频分量中阶数最低的IMF对应的能量。

10.根据权利要求9所述的一种局部放电信号的特征提取方法，其特征在于，所述的E^*的计算公式为：

其中，IMF^*为样本熵提取的高频分量中阶数最低的本征模态函数。

Images