基于广义S变换和随机子空间的局放窄带干扰抑制方法

2021-11-24宋立业蒲霄祥李希桐

电工电能新技术 2021年11期

宋立业，蒲霄祥，李希桐

(辽宁工程技术大学电气与控制工程学院，辽宁葫芦岛 125105)

1 引言

局部放电(Partial Discharge，PD)监测是有效诊断电力电缆绝缘状态的手段之一，一般情况下PD的信号能量较弱，会被现场的电磁干扰所掩盖，阻碍了PD信号的提取识别[1,2]。PD信号面临的电磁干扰主要分为随机脉冲干扰、周期性窄带干扰和白噪声三大类，其中周期性窄带干扰拥有噪声能量强和持续时间长的特点，会严重污染PD信号，因此研究PD信号中窄带干扰抑制方法具有重大的意义[3]。

目前，国内外学者针对PD信号的窄带干扰抑制方法开展了大量研究。文献[4]提出利用快速傅里叶变换配合频域中阈值抑制窄带干扰，该方法操作简单、计算速度快，但是容易受到频谱泄露的影响造成严重的边缘效应。文献[5,6]提出利用小波分解方法开展PD信号去噪，该方法可以有效分析PD信号的时频特征，进而取得较好的去噪效果，但是该方法需要人为设定小波基、分解层数和各层阈值，去噪效果受人为因素影响较大。文献[7]提出利用经验模态分解方法对PD信号进行自适应分解去噪，该方法不需要选择基函数，能更好地适应多种PD信号，但是该方法存在端点效应和模态混叠的问题。文献[8,9]提出利用奇异值分解实现窄带干扰和PD信号的分离，该方法仅需要确定奇异值阈值便可以重构出纯净的PD信号，但是该方法难以准确给出奇异值阈值，同时该方法在窄带干扰幅值较低以及干扰和PD信号的频率出现混叠时去噪效果较差。

广义S变换具有良好的时频分辨能力，能有效地分离出时频特征不同的PD信号和窄带干扰，因此被逐渐用于PD信号的窄带干扰抑制中。文献[10]提出利用广义S变换模时频矩阵的能量分布重构窄带干扰进行去噪，该方法可以有效在时频域中分离出PD信号和窄带干扰，但是没有考虑窄带干扰相位的影响。文献[11]提出利用广义S变换配合奇异值分解算法实现窄带干扰抑制，虽然解决了奇异值分解算法中奇异值阈值的选取问题，但是仍在窄带干扰和PD信号的频率出现混叠时去噪效果较差。

综上所述，本文提出一种基于广义S变换和随机子空间的局部放电窄带干扰抑制方法。该方法首先通过广义S变换得到染噪PD信号的模时频矩阵，然后根据窄带干扰和PD信号各自的时频特征确定窄带干扰数量和分离出不含PD的时间片段，接着利用随机子空间算法进行窄带干扰数据重构，最后提取出纯净的PD信号。分别利用本文方法、广义S变换模矩阵方法和频率切片小波变换方法对仿真和实测的染噪PD数据进行窄带干扰抑制，对比结果显示，相比于传统方法，本文方法能更好地抑制局部放电信号中的窄带干扰。

2 算法原理

2.1 广义S变换

S时频变换早期由Stockwell提出，该方法集合了短时傅里叶变换和连续小波变换的优点[10,11]。S变换中窗函数选用了参数可变的高斯函数，该函数的时宽和频率成反比，幅值与频率成正比，因此S变换可以在低频区间获得较好的频率分辨率，而在高频区间获得较好的时间分辨率，拥有优异的时频分析能力，对于提取PD信号这类非平稳信号的时频分布特征具有良好的效果。

信号y(t)的S变换可以表示为:

(1)

式中,f为频率；ω(t-τ,f)为高斯窗函数；t、τ为时间变量。

S变换中窗函数表达式为：

(2)

从式(2)中可以看出，S变换中窗函数的幅值和时宽是随频率变化而变化的，因此S变换中时频分辨率能够跟随频率变化而变化，能改善短时傅里叶变换中不变的时频分辨率，同时不需要考虑连续小波变换中基函数的选择问题，凭借上述优点，S变换具有良好的使用前景。

为了使S变换能够适用更多的使用场景，Pinnegar设计了调节因子λ(λ>0)对S变换中高斯窗函数进行了改进[10,11]，得到新的窗函数为：

(3)

使用式(3)作为窗函数的S变换被称为广义S变换，从式(3)中可以看出，广义S变换可以通过调节λ实现高斯窗函数和f的变化率调节。当λ<1时，高斯窗的时宽增加，幅值减少，广义S变换的频率分辨率增加，时间分辨率下降；当λ>1时，高斯窗的时宽减少，幅值增加，广义S变换的频率分辨率下降，时间分辨率增加，由此实现时频分辨率的有效调节。

由于实际信号通常为离散数据，因此需要对广义S变换做离散化处理，令f=n/(NT)和τ=iT，其中T为采样周期，N为采样个数，得到离散的广义S变换为：

(4)

式中，n，i=0，1，…，N-1；Y是y的傅里叶系数。

通过式(4)得到信号的广义S变换复数矩阵。为了能直观分析信号在时频域中的能量分布，对该复数矩阵进行求模，得到广义S变换模矩阵(Generalized S-transform Modular Matrix，GSMM)。该矩阵能直接反映信号能量的时间和频率分布特性，因此可以利用GSMM对信号的时频特征进行深度分析。

2.2 随机子空间

随机子空间算法可以有效地计算离散系统的模态参数，并且有着计算精度高、受采样条件影响小和步骤简单的优点[12]。其中基于协方差驱动的随机子空间算法凭借着优异的计算效率被广泛使用[13]，本文将其引入用于窄带干扰参数估计。

将离散信号yk以离散系统的状态方程可以表示为：

(5)

式中，A为离散系统的状态矩阵；C为离散系统的输入矩阵；xk和yk分别为离散系统k时刻的状态量和输出量；wk为离散系统中白噪声；vk为测量中白噪声。对应本文中，yk是离散的窄带干扰信号。

将yk建立三个时间矩阵分别为：

(6)

(7)

(8)

式中，Yp为“过去”时间矩阵；Yf1为第一个“未来”时间矩阵；Yf2为第二个“未来”时间矩阵；b为时间矩阵的列数，该值越大，说明离散系统输出数据量越多，此时各参数估计越精准，但是yk的数据量是有限的，因此本文取b>10(N-b-1)，得到b为：

b=Floor(10(N-1)/11)

(9)

式中，Floor是向下取整。

将三个时间矩阵Yp、Yf1和Yf2构建托普利茨矩阵T1和T2为：

T1=Yf1YpT

(10)

T2=Yf2YpT

(11)

接着将T1开展奇异值分解为：

(12)

式中，S1、S2分别为信号和噪声主导的奇异值对角矩阵；U1、U2分别为信号和噪声主导的左单位奇异矩阵；V1、V2分别为信号和噪声主导的右单位奇异矩阵。

从统计理论上分析，wk和vk应该是无相关性的两组噪声，由此根据离散系统的理论可得：

(13)

式中，Oa为观测矩阵；G=E(x(k+1)ykT)，E为期望；Гa为控制矩阵；a=N-b-1。

将式(12)和式(13)进行联合分析可得：

(14)

再依据式(13)可得：

T2=OaAΓa

(15)

进一步联合式(14)和式(15)得到A为：

(16)

对A开展特征值分解为：

A=ψΛψ-1

(17)

式中，Λ=diag(zs)，s=1，2，…，p；zs是A的第s个特征值；p为式(5)中系统的阶数，对应本文为窄带干扰个数的2倍；ψ为A的特征矩阵。

利用最小二乘法求解下式即可得到参数Bs。

(18)

利用zs和Bs得到窄带干扰的幅值Qs、频率fs和相位θs分别为：

Qs=2|Bs|

(19)

(20)

(21)

由于窄带干扰为实信号，因此确定fs>0的对应数据组为真实数据，进而估计得到窄带干扰的幅值、频率和相位。

2.3 PD窄带抑制流程

本文结合广义S变换和随机子空间提出PD窄带干扰抑制方法的具体步骤为：

(1)将染噪PD信号开展广义S变换，得到对应的广义S变换模矩阵GSMM。

(2)借助窄带干扰时间分布长，频率能量集中的特点确定窄带干扰个数。

(3)借助PD信号时间分布短，频率能量分布宽的特点确定染噪PD信号中无PD的最长时间片段。

(4)利用随机子空间算法对染噪PD信号中无PD的最长时间片段进行处理，估计窄带干扰参数。

(5)利用估计的窄带干扰参数值对窄带干扰信号进行重构，在时域中将染噪PD信号减去窄带干扰的重构值，得到窄带干扰抑制后的PD信号。

3 仿真测试

3.1 仿真局放模型

PD信号通常呈现衰减振荡的特性，因此可以用如式(22)和式(23)所示的单指数衰减振荡型和双指数衰减振荡型信号模型进行模拟[14]。

g1(t)=Ce-t/ηsin(2πfgt)

(22)

g2(t)=C(e-1.3t/η-e-2.2t/η)sin(2πfgt)

(23)

式中，C为脉冲幅值；η为衰减系数；fg为振荡频率。

本节模拟4个PD脉冲信号，各脉冲信号参数如表1所示，其中脉冲Ⅰ和Ⅲ是单指数衰减振荡型脉冲，脉冲Ⅱ和Ⅳ是双指数衰减振荡型脉冲，采样频率为30 MHz。

表1 仿真PD脉冲参数

仿真中窄带干扰由4个不同频率、幅值和初始相位的正弦波叠加组成，根据文献[11、15]的仿真窄带干扰参数范围，本文将仿真中窄带干扰的参数值设置如表2所示。由此得到纯净的PD信号、窄带干扰信号和染噪PD信号如图1所示。从图1中可以看出，在染噪PD信号中，由于窄带干扰的存在，PD信号几乎完全被淹没，难以直接进行识别和提取。

表2 仿真窄带干扰参数

图1 仿真的PD信号

3.2 局放去噪

由于广义S变换仍属于短时傅里叶变换，因此其时间分辨率和频率分辨率并不能同时达到最高，本文将λ取为0.4以同时获得较好的时间分辨率和频率分辨率。通过广义S变换处理图1(c)中的染噪PD信号，得到染噪PD信号的GSMM如图2所示。从图2中可以看出窄带干扰信号和PD信号存在明显不同的特征：窄带干扰信号时间上分布较长，频率上分布较集中；PD信号时间上分布较短，频率上分布较宽。

以此可以通过分析图2中染噪PD信号的GSMM确定窄带干扰的个数为4，同时可以在时间横轴上进行区域划分，得到无PD区域。对于图2而言，横轴上采样点数为1 400～2 400区域可以被视为无PD区域，该区域中主导信号为窄带干扰信号。将最长时段的无PD区域作为窄带干扰参数估计时间段，选择图1(c)中染噪PD信号中该时间段的数据进行随机子空间算法处理，得到窄带干扰参数估计值如表3所示。对比表2和表3可以看出，该方法可以有效地估计窄带干扰的各参数值。值得说明的是，由于仅截取了染噪PD信号中部分时段区段进行参数估计，因此估计得到的相位和原始参数的相位是不一致的。

表3 仿真染噪PD信号的窄带干扰参数估计值

利用表3中窄带干扰参数估计值重构整个时间段的窄带干扰，重构的窄带干扰Dg如下：

(24)

式中，t1为截取区段的起始时刻。

将染噪PD信号减去重构的窄带干扰得到去噪结果如图3(a)所示，为了进行对比，利用广义S变换模矩阵方法[10]和频率切片小波变换方法[16]处理图1(c)中染噪PD信号，得到去噪结果分别如图3(b)和图3(c)所示。从图3的对比结果可以看出，对于广义S变换模矩阵方法而言，该方法没有考虑窄带干扰相位的影响，因此对于本文中窄带干扰相位不为0的情况，该方法无法实现窄带干扰抑制；对于频率切片小波变换方法而言，去噪后波形存在较大残余噪声，同时波形中存在明显的边缘效应，整体去噪效果较差；本文方法去噪后波形的噪声较小，PD信号的细节得到保留，利于后续PD信号的提取识别分析。

图3 仿真染噪PD信号的窄带干扰抑制结果

为了进一步说明3种方法的去噪效果，本文引入信噪比SNR、均方误差MSE和波形相似度NCC三个评价参数开展分析[17]，其具体计算公式为：

(25)

(26)

(27)

式中，g(n)为原始PD信号；d(n)为各方法去噪后的信号。其中SNR越大，MSE越小，NNC越接近于1时，去噪效果越好。

计算得到图3中各方法去噪后波形的SNR、MSE和NCC如表4所示，从表4中可以看出，和广义S变换模矩阵方法、频率切片小波变换方法相比，本文方法能更好地去除染噪PD信号中窄带干扰，去噪后波形更好地保留了原始PD信号的波形特征，利于后续的分析研究。

表4 窄带干扰抑制结果对比

4 实测研究

为了验证本文方法在实际测试中的可行性，在实验室中对10 kV电缆开展工频局放测试，电缆终端头中制作有半导电层突刺缺陷[18]，监测方法为高频电流法[7]，采样率设置为200 MHz，其具体接线图如图4所示。

图4 工频局放测试接线图

由于实验室中采集的PD信号噪声较小，因此在采集的PD信号中添加多个窄带干扰信号。文献[19]指出高频电流法受到的窄带干扰主要包括中波段0.5～1.6 MHz、短波段 2.3～25 MHz和调频段88～108 MHz的广播信号。因此本文在此基础上随机选择窄带干扰信号的幅值分别为2 mV、4 mV、3 mV和2 mV；频率分别为5.06 MHz、10 MHz、15.18 MHz和22 MHz；相位分别为π/4 rad、π/3 rad、π/5 rad和π/2 rad，得到实测的染噪PD信号如图5所示，从图5中可以看出，由于窄带干扰的污染，PD信号无法直接被识别提取。

图5 实测的染噪PD信号

利用本文方法对图5中染噪PD信号进行去噪处理，得到染噪PD信号的GSMM图如图6所示，借助PD信号和窄带干扰的时频特征，可以确定图6中窄带干扰数目为4，同时很容易在时间横轴上确定无PD区域。对于图6而言，时间横轴上采样点数为600～1 400区域可以被视为无PD区域，将图5中染噪PD信号对应时间段的数据利用随机子空间算法进行窄带干扰参数估计，得到窄带干扰参数估计结果如表5所示，利用该结果对窄带干扰进行重构，得到去噪后波形如图7(a)所示。为了进一步说明本文方法的优越性，再利用广义S变换模矩阵方法和频率切片小波变换方法对图5中染噪PD信号进行去噪，得到对应去噪结果如图7(b)和图7(c)所示。从图7的对比结果中可以看出，广义S变换模矩阵方法抑制窄带干扰失败，没有提取出局放波形；频率切片小波方法的去噪结果中存在明显的窄带干扰抑制不干净现象，同时存在明显的边缘效应，因此上述2种传统方法的去噪效果较差。相对于传统方法而言，本文方法可以更加有效地提取出PD波形，残余噪声更小。