结合广义S变换和快速独立分量分析的局放信号中窄带干扰抑制方法

2022-10-17李向群伍亚萍庄旭菲吉飞敏张燕

现代电力 2022年5期

李向群，伍亚萍，庄旭菲，吉飞敏，张燕

（1. 西北民族大学电气工程学院, 甘肃省兰州市 730030；2. 内蒙古工业大学信息工程学院,内蒙古自治区呼和浩特市 010051）

0 引言

电力电缆是输配电系统的重要组成部分[1]，而局部放电（partial discharge，PD）是电力电缆绝缘劣化的早期表现形式之一，因此PD信号的监测和分析已经成为电力电缆绝缘评估的有效手段[2]。但是，由于实际采集的PD信号能量较小，而实际现场的电磁干扰能量较大，所以PD信号通常会被电磁干扰淹没，严重妨碍PD信号的监测和分析。常见的电磁干扰主要为随机脉冲干扰[3]、白噪声[4-5]和周期性窄带干扰[6-7]。其中周期性窄带干扰有着较强的能量和较长的持续时间，通常会将采集的PD信号完全淹没，因此PD信号中周期性窄带干扰的抑制技术具有重大的研究意义[8]。

在PD信号中周期性窄带干扰的抑制技术上，国内外学者已经开展了大量的研究工作。文献[8-9]利用快速傅里叶变换（fast Fourier transform，FFT）方法有效削弱了周期性窄带干扰的能量，但是由于频谱泄露问题的存在，所以该类方法难以获得较好的干扰抑制效果。文献[10]利用小波分解方法的优异时频分析能力提升了干扰抑制效果，但是由于PD信号波形存在多样性的特点，所以该类方法难以选择合适的小波基函数。文献[11]利用经验模态分解方法的自适应分解能力对窄带干扰进行抑制，避免了基函数的选取问题，但是由于端点效应和模态混叠问题的存在，所以该类方法的运用范围极其有限。文献[12-13]利用奇异值分解方法将信号分为窄带干扰子空间和PD信号子空间，该类方法利用奇异值分解方法的信号自适应分解能力实现了窄带干扰的抑制，但是该类方法存在奇异值阈值选取困难和部分窄带干扰抑制效果差的问题。

为了解决现有技术存在的问题，本文提出基于广义S变换和快速独立分量分析的局放信号窄带干扰抑制方法。该方法首先利用广义S变换获取窄带干扰的特征区域，然后在窄带干扰的特征区域中，利用Candan算法对窄带干扰频率进行精准估计，最后利用改进的快速独立分量分析方法分离PD信号，从而抑制PD信号中窄带干扰。仿真和实际的测试结果说明，相比于传统方法，本文方法能够有效提升窄带干扰的抑制效果。

1 窄带干扰抑制方法

1.1 仿真信号

大量的研究结果表明，PD信号可以采用单指数振荡衰减函数和双指数振荡衰减函数进行模拟分析研究[14]，即仿真PD信号的函数表达式为

式中：fp是PD信号的振荡频率；τ是PD信号的衰减系数；A是PD信号的幅值系数。

通过式（1）和式（2）的函数模型和表1中仿真PD脉冲的参数值，得到仿真PD脉冲如图1（a）所示，其中脉冲1和脉冲2为单指数振荡衰减函数模型，脉冲3和脉冲4为双指数振荡衰减函数模型。该仿真中采样时长设置为50 μs，采样频率fs设置为50 MHz。

表1 仿真PD脉冲的参数Table 1 Parameters of simulated PD pulse

图1 仿真的PD信号Fig. 1 Simulated PD signal

一般情况下，周期性窄带干扰会表现为正弦信号或余弦信号，因此通常可以采用三角周期函数进行模拟[15]，具体的函数表达式为

式中：m是窄带干扰的总数；Bi是窄带干扰的幅值；fni是窄带干扰的频率；θi是窄带干扰的相位。

本文将m设置为3，并且各仿真窄带干扰的设置值如表2所示，以此得到仿真的窄带干扰波形如图1（b）所示。并将该窄带干扰叠加到仿真PD脉冲中，得到染噪PD信号波形如图1（c）所示，从图1（c）中可以看出，原始的PD信号已经完全被窄带干扰所掩盖，无法进行识别分析，因此需要进行窄带干扰抑制。

表2 仿真窄带干扰的参数Table 2 Parameters of simulated narrowband interference

1.2 广义S变换

S变换是由Stockwell研究得到的一种信号时频分析方法，该方法同时兼顾了短时傅里叶变换和小波变换的优点[16-17]。广义S变换在S变换的基础上引入了调节因子λ，以此实现时频分辨率的人为调节，获得更加适宜的时频分辨率[18]。

定义广义S变换的窗函数w(η-t,f,λ)为

从式（4）中可以看出，通过人为设定λ值可以改变窗口长度，当λ>1时，可以减小窗口长度，从而提升时域分辨率；当λ<1时，可以增加窗口长度，从而提升频率分辨率；当λ=1时，广义S变换退化为S变换。

对于信号x的离散形式x(i)，i=0，1，…，N-1而言，得到对应的广义S变换为

式中：T是采样周期；n=0，1，…，N-1；N是采样总长度。

通过式（5）可以得到x(i)的广义S变换时频分析矩阵，该矩阵中行向量代表时间；列向量代表频率；矩阵中元素为复数，包含了广义S变换结果的幅值信息和相位信息。为了方便对广义S变换结果进行分析，对广义S变换时频分析矩阵进行求模处理得到模矩阵为

通过对式（6）中模矩阵进行分析，可以获取信号中各频率分量的时间分布情况。

1.3 窄带干扰特征量分析

通过广义S变换可以对染噪PD信号进行时频分析，从而在时频分布图谱上对窄带干扰和PD信号进行单独分析。由于S变换的时频分辨率完全由分析频率所决定，无法进行人为调节，所以将广义S变换的λ设置为0.4，以获得合适的时间分辨率和频率分辨率，便于在时频分布图谱中分析窄带干扰和PD信号的时频特征，从而提取出窄带干扰的特征区域。

图2为图1（c）中染噪PD信号的广义S变换模矩阵，即染噪PD信号的时频分布图谱。通过对染噪PD信号的时频分布图谱进行分析，可以划分出窄带干扰的特征区域，从而准确提取窄带干扰的特征量，实现窄带干扰抑制。在染噪PD信号的时频分布图谱中，周期性窄带干扰的时间分布较长，频率分布较集中；PD信号的时间分布较短，频率分布较广，因此可以很容易地在时频分布图谱中区分出窄带干扰和PD信号。在区分出窄带干扰和PD信号后，可以确定窄带干扰的总数为3，然后对各窄带干扰对应的时频矩阵区域进行分析，将没有和PD信号发生混叠的时频子矩阵区域划出，作为各窄带干扰的特征区域，以此得到窄带干扰特征区域Ⅰ、Ⅱ和Ⅲ。

图2 仿真PD信号的时频分布Fig. 2 Time frequency distribution picture of simulated PD signal

为了在各窄带干扰特征区域中准确估计窄带干扰的频率，本文引入了Candan算法[19]。该算法在傅里叶频率估计算法的基础上，考虑了原有算法存在的栅栏效应问题，能取得较高的频率估计精度，同时该算法避免了加窗赋权和多次迭代的计算步骤，具有较低的计算复杂度，利于工程中实现，其具体计算步骤如下：

1）对窄带干扰特征区域进行边缘标记，得到窄带干扰特征区域的频率范围为[fmin，fmax]，采样点范围为[nmin，nmax]，得到采样点范围内长度N1=nmax-nmin+1，同时获得离散频率点范围为[Floor(N1·fmin/fs)，Ceil(N1·fmax/fs)]，Floor(*)是指向下取整，Ceil(*)是指向上取整，以此得到窄带干扰的特征采样点内数据为y(i)，i=0，1，…，N1-1，并且y(i)=x(i+nmin)。

2）对y(i)中离散频率点范围[Floor(N1·fmin/fs)，Ceil(N1·fmax/fs)]进行离散傅里叶变换，得到傅里叶系数为

3）确定Y(k)的幅值最大值位置，并记录对应的离散频率点为p，对应的傅里叶系数为Y(p)，同时确定附近的傅里叶系数分别为Y(p-1)和Y(p+1)。

4）计算频率修正值为

5）确定窄带干扰的精准频率为

利用上述方法分别对窄带干扰特征区域Ⅰ、Ⅱ和Ⅲ中窄带干扰进行频率估计，得到窄带干扰频率估计值如表3所示。从表3中可以看出，本文所提Candan算法可以精准估计各窄带干扰的频率，说明本文方法可以有效解决傅里叶频率估计算法中栅栏效应的问题，同时该方法无需加窗处理和多次迭代计算，执行效率较高。

表3 仿真PD信号的窄带干扰频率估计结果Table 3 Narrowband interference frequency estimation of simulated PD signal

1.4 快速独立分量分析

染噪PD信号可以被视为PD信号和窄带干扰的线性混合，进一步可以被视为PD信号、正弦信号和余弦信号的线性混合，因此本文采用快速独立分量分析算法对染噪PD信号进行分离，该算法具体原理可参考文献[20]。

由于快速独立分量分析算法在开展混合信号的分离时，必须要求混合信号的数目大于源信号的数目，因此为了保证快速独立分量分析算法可以有效分离出原始PD信号，需要利用1.3节频率估计结果对混合信号矩阵进行扩展。对于每个窄带干扰，构建2个混合信号分别为sin(2π)和cos(2π)，和染噪PD信号共同构成混合信号矩阵。

由于快速独立分量分析算法自身存在的缺陷，所以分离的信号存在顺序、幅值和极性的随机性，即在分离的信号中难以确定哪组信号为PD信号，同时分离的PD信号和原始PD信号的幅值和极性可能存在差异。为了解决上述问题，本文在染噪PD信号前端添加长度为300的校正信号z(i)，该校正信号为

式中max(*)是求数组的最大值。

图3 混合信号Fig. 3 Mixed signal

对于分离信号中PD信号的判断问题，本文引入绝对相关数作为评判标准，定义各组分离信号前端长度为300的信号为r(i)，得到r(i)和z(i)的绝对相关数ρ为

取绝对相关数最大的一组分离信号作为分离后的PD信号(i)，然后利用式（11）对信号(i)进行修正，使得窄带干扰抑制后的PD信号q(i)更接近于原始PD信号。

2 PD信号的窄带干扰抑制

2.1 干扰抑制流程

本文所提PD信号中窄带干扰抑制方法的具体步骤为：

1）对染噪PD信号进行广义S变换，利用PD信号和窄带干扰的时频分布差异，在广义S变换中划分出窄带干扰特征区域。

2）在各窄带干扰的特征区域中，利用Candan算法对窄带干扰频率进行精准估计。

3）利用估计得到的窄带干扰频率，构建混合信号sin(2π)和cos(2π)。

4）在染噪PD信号的前端添加校正信号，并在混合信号sin(2π)和cos(2π)的前端添加零值信号，并将上述信号共同构成混合信号矩阵。

5）在利用快速独立分量分析算法对混合信号矩阵进行分析后，根据绝对相关数确定哪组信号为PD信号，同时利用校正信号修正分离的PD信号，得到窄带干扰抑制结果。

2.2 干扰抑制效果

利用本文所提PD信号中窄带干扰抑制方法对图1（c）中染噪PD信号进行去噪，通过图2的时频分布图可以确定窄带干扰的数目为3，对应的窄带干扰特征矩阵也为3个，以此得到窄带干扰抑制效果如图4（a）所示。

为了说明本文方法的优越性，使用奇异值分解方法[21]和FFT滤波方法[10]对图1（c）中染噪PD信号进行窄带干扰抑制，得到干扰抑制结果分别如图4（b）和图4（c）所示。从图4中可以看出，本文所提方法可以有效抑制窄带干扰，并且窄带干扰抑制后的PD信号波形和原始PD信号波形的幅值和极性保持一致，说明本文方法可以有效恢复PD信号的波形特征，利于PD信号的后续分析。而传统的奇异值分解方法和FFT滤波方法的窄带干扰抑制结果中存在明显的残余噪声，窄带干扰抑制效果较差。

为了量化说明图4中各窄带干扰抑制方法的效果，本文引入3组参数对去噪效果进行评价，分别为信噪比系数SNR、波形相似系数NCC和均方误差系数MSE[13]：

图4 仿真PD信号的窄带干扰抑制结果Fig. 4 Suppression results of narrowband interference for PD signal

计算得到图4中各窄带干扰抑制结果的评价参数值如表4所示。对比表4中各方法的评价参数可以看出，在3种方法中，本文所提方法能取得最高的SNR和NCC，并且能取得最低的MSE，说明本文所提方法可以很好恢复得到PD信号的波形，波形的畸变失真较少，利于后续的PD信号分析。

表4 窄带干扰的抑制效果Table 4 Suppression effect of narrowband interference

需要说明的是，本文方法的计算复杂度较高，但是窄带干扰的抑制效果更好，因此本文方法更加适合PD信号的离线分析。

3 实测信号的窄带干扰抑制

为了分析本文所提方法的实际运用效果，在实验室中制作缺陷电缆样本，并开展工频局部放电测试。电缆样本选用10 kV交联聚乙烯电缆，设置的缺陷为终端头刀痕缺陷，具体的试验接线图如图5所示。图5中示波器的采样率设置为100 MHz，采样点数为1300。由于实验室中电磁干扰较小，因此向采集得到的PD信号中叠加频率分别为1.34、3.38和10.43 MHz，幅值分别为3、2和3 mV，相位分别为π/2、π/6和π/3 rad的窄带干扰信号，得到实测的染噪PD信号如图6所示。从图6中可以看出，在窄带干扰的影响下，PD信号已经无法被识别提取，所以窄带干扰的抑制是必要的。

图5 局放测试接线示意图Fig. 5 Wiring diagram of PD measurement

利用广义S变换分析图6中实测的染噪PD信号，得到染噪PD信号的时频分析图谱如图7所示，从图7中可以看出，窄带干扰的数目为3，同时划定各窄带干扰对应的时频区域，得到窄带干扰特征区域Ⅰ、Ⅱ和Ⅲ。在窄带干扰特征区域内，对窄带干扰的频率进行估计，得到实测染噪PD信号中各窄带干扰的频率如表5所示。

从表5中可以看出，本文所提方法可以精准估计实测染噪PD信号的窄带干扰频率，利于后续窄带干扰的抑制，继续使用本文方法对窄带干扰进行抑制，得到窄带干扰抑制结果如图8（a）所示。同时为了对比传统的方法，进一步采用奇异值分解法和FFT滤波法对图6中实测染噪PD信号进行窄带干扰抑制，分别得到干扰抑制结果如图8（b）和图8（c）所示。对比图8中各方法的窄带干扰抑制结果可以看出：奇异值分解方法和FFT滤波方法的窄带干扰抑制结果中存在明显的残余噪声（图8中圈内标记所示），窄带干扰抑制效果有限；本文所提方法的窄带干扰抑制效果明显优于传统方法，能很好地抑制窄带干扰。