基于时频谱图分离的局部放电窄带干扰抑制方法

摘要

本发明公开了一种基于时频谱图分离的局部放电窄带干扰抑制方法，包括：采用广义S变换获取染噪PD信号的二维时频复数矩阵，然后求模得到染噪PD信号时频谱图；利用数学形态学分离出染噪PD信号时频谱图中窄带干扰和原始PD信号的时频谱图，以提取窄带干扰和PD信号的特征量；将染噪PD信号划定为信号帧和仅含噪声的噪声帧，采用总体最小二乘不变旋转矢量技术分析噪声帧内时域数据，以提取窄带干扰的时域特征量；重构出窄带干扰的时域波形，从染噪PD信号的时域波形中去除窄带干扰的时域波形。本发明应用时能提升窄带干扰影响下PD信号时频分析的准确率，从而有效抑制PD信号中窄带干扰。

权利要求

1.基于时频谱图分离的局部放电窄带干扰抑制方法，其特征在于，包括：S1、采样获取染噪PD信号；

S2、采用广义S变换对染噪PD信号的离散信号进行时频变换处理，得到二维时频复数矩阵，其中，二维时频复数矩阵中横轴为时间采样点，纵轴为频率采样点；

S3、对二维时频复数矩阵求模得到关于染噪PD信号的时频谱图；

S4、采用数学形态学获取染噪PD信号时频谱图中各行分量的平滑分量和非平滑分量，将分离出的平滑分量和非平滑分量作为行分量构成新的矩阵，分离出窄带干扰和原始PD信号的时频谱图，以确定窄带干扰的数目与PD脉冲的时间区域；

S5、将采样得到的染噪PD信号根据时间轴划定为含有PD信号的信号帧和仅含噪声的噪声帧，采用总体最小二乘不变旋转矢量技术分析噪声帧内时域数据，以提取窄带干扰的时域特征量；

S6、重构出窄带干扰的时域波形，从染噪PD信号的时域波形中去除窄带干扰的时域波形，得到降噪后PD信号；

所述步骤S2中采用广义S变换对染噪PD信号的离散信号进行时频变换处理包括以下步骤：染噪PD信号x(t)的广义S变换结果G(τ,f,λ)被定义为式中，t和τ为两组时间变量；f是频率；λ是调节因子；w(t‑τ,f,λ)是高斯窗函数，对应的表达式为：令x(n)为染噪PD信号x的离散信号，同时定义f＝n/(NT)、τ＝iT，其中T是采样周期，N是x(n)的数据总量，得到x(n)的广义S变换结果为：式中，i、m、n是x(n)的广义S变换结果中三组变量，被定义为0，1，…，N‑1；

利用x(n)的广义S变换结果对x(n)开展时频变换处理后得到二维的时频复数矩阵；

所述步骤S4中采用数学形态学获取染噪PD信号时频谱图中各行分量的平滑分量和非平滑分量包括以下步骤：所述染噪PD信号时频谱图中各行分量Gt(nGt)的膨胀变换结果公式为：(Gt⊕g)(nGt)＝

max{Gt(nGt‑ig)+g(ig)|(nGt‑ig)∈DGt且ig∈Dg}腐蚀变换结果公式为：

(Gt⊙g)(nGt)＝

min{Gt(nGt+ig)‑g(ig)|(nGt+ig)∈DGt且ig∈Dg}式中，g(ig)是结构元素；⊕是膨胀运算符；⊙是腐蚀运算符；DGt是Gt(nGt)的定义域；Dg是g(ig)的定义域，nGt是数组Gt中某个索引，ig是数组g中某个索引；

通过膨胀变换结果公式和腐蚀变换结果公式的级联组合，得到数学形态学中开运算和闭运算分别为：(Gt·g)(nGt)＝((Gt⊕g)⊙g)(nGt)式中，是开运算；·是闭运算；

将开、闭运算进行级联组合，分别得到形态开‑闭滤波器运行结果[Foc(Gt)](nGt)和形态闭‑开滤波器运行结果[Fco(Gt)](nGt)为:采用形态开‑闭滤波器和形态闭‑开滤波器的混合运算，得到信号Gt(nGt)中的平滑分量[F1(Gt)](nGt)和非平滑分量[F2(Gt)](nGt)分别为：[F1(Gt)](nGt)＝([Foc(Gt)](nGt)+[Fco(Gt)](nGt))/2[F2(Gt)](nGt)＝Gt(nGt)‑[F1(Gt)](nGt)；

所述步骤S5中采用总体最小二乘不变旋转矢量技术分析噪声帧内时域数据包括以下步骤：S51、将噪声帧内时域数据信号y(ny)构成Hankle矩阵H为：式中，Ny是信号y(ny)的数据长度，L是H的列数；

S52、将H开展奇异值分解，得到

T

H＝UΣV

T

式中，是取矩阵的共轭转置；U和V分别是左、右正交矩阵；Σ是对角矩阵，其中的对角元素是矩阵H的奇异值；

S53、分离出矩阵U中第1行到第L行，第1列到第P列的子矩阵为矩阵U1；同时分离出矩阵U中第2行到第L+1行，第1列到第P列的子矩阵为U2，以此得到U1和U2分别为：U1＝U[1:L,1:P]

U2＝U[2:L+1,1:P]

式中，P是分离后时频图谱确定的窄带干扰个数的2倍；

S54、构建矩阵Z＝[U1 U2]，然后对Z开展奇异值分解，得到式中： 和 均是矩阵Z奇异值分解后的单位特征矩阵； 是矩阵Z奇异值分解后的对角矩阵；V11、V12、V21和V22是 的子矩阵；

S55、利用V12和V22构建矩阵ψ为‑1

ψ＝‑V12V22

S56、对ψ开展特征值分解，得到对应的特征值为μk(k＝1,2,…,P)，从而估计得到各窄带干扰分量的频率 为：式中，arg(*)是计算复数的相角；

S57、利用最小二乘法计算各窄带干扰分量的幅值 和相角 具体计算公式为：Y＝μc

T

Y＝[y(0)，y(1)，…，y(Ny‑1)]式中，c是中间变量矩阵，ck是c中元素。

2.根据权利要求1所述的基于时频谱图分离的局部放电窄带干扰抑制方法，其特征在于，所述信号帧内窄带干扰相位 通过相角 开展移相处理得到，其计算公式为：式中，ΔT是信号帧和噪声帧的时延。

说明书

基于时频谱图分离的局部放电窄带干扰抑制方法

技术领域

[0001] 本发明涉及高压电力设备绝缘状态诊断技术，具体是基于时频谱图分离的局部放电窄带干扰抑制方法。

背景技术

[0002] 高压电力设备的绝缘状态可通过监测局部放电(Partial Discharge,PD)信号进行诊断，通常现场采集的PD信号会受到电磁干扰影响，导致难以直接开展PD信号的分析及特征提取。PD信号面临的电磁干扰主要为脉冲型干扰、白噪声干扰及周期性窄带干扰，其中，周期性窄带干扰拥有持续时间长、能量强和随机性大的显著特点，导致PD信号容易被周期性窄带干扰完全淹没。因此，为了利用PD信号对高压电力设备的绝缘状态诊断时提升精确性，对周期性窄带干扰的抑制尤其重要。

[0003] 目前，国内外学者对PD信号中周期性窄带干扰的抑制技术开展了大量研究。樊高辉、刘尚合、刘卫东及王雷于2017年04月06日在《高电压技术》发表的名称为“FFT谱最小熵解卷积滤波抑制放电信号中的周期性窄带干扰”的文献，其在快速傅里叶变换(fast Fourier transform，FFT)算法的基础上，引入了解卷积滤波法和经典阈值法用于提取窄带干扰的频点，然后将窄带干扰对应频点的FFT系数进行压缩，达到削弱窄带干扰能量的目的。虽然该方法可以有效削弱窄带干扰的频域能量，但是由于FFT算法存在频谱泄露等缺陷和阈值法存在阈值选取困难的问题，所以该方法最终的降噪效果不够理想。马星河和张登奎于2021年08月17日在《电工技术学报》发表的名称为“基于改进经验小波变换的高压电缆局部放电噪声抑制研究”的文献，其借助小波分解算法的优异时频分析能力，能一定程度上削弱窄带干扰的能量，但是该类方法难以选择适宜的小波基函数、小波分解层数和小波分解阈值。魏海增、马宏忠、黄涛及黄烜城于2019年05朋15日在《电力系统及其自动化学报》发表的名称为“基于EMD的ICA降噪方法在电厂开关柜局部放电信号中的应用”的文献，其使用经验模态分解算法对PD信号中窄带干扰进行分离，但是该方法的稳定性较差，容易出现端点效应、模态混叠等问题。徐永干、姜杰、唐昆明、张太勤、罗建及谢敏于2019年09月17日在《电网技术》发表的名称为“基于Hankel矩阵和奇异值分解的局部放电窄带干扰抑制方法”的文献，毕潇文、钟俊、张大堃、周电波及阮莹于2021年03月25日在《电网技术》发表的名称为“基于改进奇异值与经验小波分解的局放去噪算法”的文献，以及杨晓丽、黄宏光、舒勤、张大堃及周电波于2020年11月25日在《高电压技术》发表的名称为“基于SVD和低秩RBF神经网络的局部放电信号提取方法”的文献，这三篇文献均提出了自适应奇异值分解法抑制染噪PD信号中窄带干扰，但是该方法存在奇异值阈值选取困难的问题，难以有效抑制小幅值的窄带干扰。宋立业、蒲霄祥及李希桐于2021年11月23日在《电工电能新技术》发表的名称为“基于广义S变换和随机子空间的局放窄带干扰抑制方法”的文献，其利用广义S变换算法得到染噪PD信号的时频谱图，然后直接在该时频谱图中确定窄带干扰和PD信号对应的时频区域子矩阵，以提取窄带干扰的特征参数，实现窄带干扰抑制，但是在窄带干扰强能量特点的影响下，PD信号的时频特征可能会被窄带干扰所掩盖，导致该类方法的可靠性下降。

发明内容

[0004] 本发明的目的在于解决现有技术抑制染噪PD信号中窄带干扰的不足，提供了一种基于时频谱图分离的局部放电窄带干扰抑制方法，其能分离出染噪PD信号时频谱图中窄带干扰和PD信号，提升窄带干扰影响下PD信号时频分析的准确率，从而有效抑制染噪PD信号中窄带干扰。

[0005] 本发明的目的主要通过以下技术方案实现：

[0006] 基于时频谱图分离的局部放电窄带干扰抑制方法，包括：

[0007] S1、采样获取染噪PD信号；

[0008] S2、采用广义S变换对染噪PD信号的离散信号进行时频变换处理，得到二维时频复数矩阵，其中，二维时频复数矩阵中横轴为时间采样点，纵轴为频率采样点；

[0009] S3、对二维时频复数矩阵求模得到关于染噪PD信号的时频谱图；

[0010] S4、采用数学形态学获取染噪PD信号时频谱图中各行分量的平滑分量和非平滑分量，将分离出的平滑分量和非平滑分量作为行分量构成新的矩阵，分离出窄带干扰和原始PD信号的时频谱图，以确定窄带干扰的数目与PD脉冲的时间区域；

[0011] S5、将采样得到的染噪PD信号根据时间轴划定为含有PD信号的信号帧和仅含噪声的噪声帧，采用总体最小二乘不变旋转矢量技术分析噪声帧内时域数据，以提取窄带干扰的时域特征量；

[0012] S6、重构出窄带干扰的时域波形，从染噪PD信号的时域波形中去除窄带干扰的时域波形，得到降噪后PD信号。

[0013] 本发明针对高压电力设备局部放电信号受周期性窄带干扰影响大的问题，提出了一种基于时频谱图分离的局部放电窄带干扰抑制方法，其结合广义S变换和数学形态学，从染噪PD信号中分离出原始PD信号和窄带干扰信号的时频谱图，从而在时频域内准确辨识出窄带干扰和原始PD信号，避免强能量的窄带干扰掩盖原始PD信号时频特征的问题，提升窄带干扰影响下PD信号时频分析的准确率。本发明再利用总体最小二乘不变旋转矢量技术准确估计窄带干扰的特征参数，以抑制染噪PD信号中窄带干扰。

[0014] 进一步的，所述步骤S2中采用广义S变换对染噪PD信号的离散信号进行时频变换处理包括以下步骤：

[0015] 染噪PD信号x(t)的广义S变换结果G(τ,f,λ)被定义为

[0016]

[0017] 式中，t和τ为两组时间变量；f是频率；λ是调节因子；w(t‑τ,f,λ)是高斯窗函数，对应的表达式为：

[0018]

[0019] 令x(n)为染噪PD信号x的离散信号，同时定义f＝n/(NT)、τ＝iT，其中T是采样周期，N是x(n)的数据总量，得到x(n)的广义S变换结果为：

[0020]

[0021] 式中，i、m、n是x(n)的广义S变换结果中三组变量，被定义为0，1，…，N‑1；

[0022] 利用x(n)的广义S变换结果对x(n)开展时频变换处理后得到二维的时频复数矩阵。

[0023] 进一步的，所述步骤S4中采用数学形态学获取染噪PD信号时频谱图中各行分量的平滑分量和非平滑分量包括以下步骤：

[0024] 所述染噪PD信号时频谱图中各行分量Gt(nGt)的膨胀变换结果公式为：

[0025]

[0026] 腐蚀变换结果公式为：

[0027]

[0028] 式中，g(ig)是结构元素；⊕是膨胀运算符；⊙是腐蚀运算符；DGt是Gt(nGt)的定义域；Dg是g(ig)的定义域，nGt是数组Gt中某个索引，ig是数组g中某个索引；

[0029] 通过膨胀变换结果公式和腐蚀变换结果公式的级联组合，得到数学形态学中开运算和闭运算分别为：

[0030]

[0031] (Gt·g)(nGt)＝((Gt⊕g)⊙g)(nGt)

[0032] 式中，是开运算；·是闭运算；

[0033] 将开、闭运算进行级联组合，分别得到形态开‑闭滤波器运行结果[Foc(Gt)](nGt)和形态闭‑开滤波器运行结果[Fco(Gt)](nGt)为:

[0034]

[0035]

[0036] 采用形态开‑闭滤波器和形态闭‑开滤波器的混合运算，得到信号Gt(nGt)中的平滑分量[F1(Gt)](nGt)和非平滑分量[F2(Gt)](nGt)分别为：

[0037] [F1(Gt)](nGt)＝([Foc(Gt)](nGt)+[Fco(Gt)](nGt))/2

[0038] [F2(Gt)](nGt)＝Gt(nGt)‑[F1(Gt)](nGt)。

[0039] 进一步的，所述步骤S5中采用总体最小二乘不变旋转矢量技术分析噪声帧内时域数据包括以下步骤：

[0040] S51、将噪声帧内时域数据信号y(ny)构成Hankle矩阵H为：

[0041]

[0042] 式中，Ny是信号y(ny)的数据长度，L是H的列数；

[0043] S52、将H开展奇异值分解，得到

[0044] H＝UΣVT

[0045] 式中，T是取矩阵的共轭转置；U和V分别是左、右正交矩阵；Σ是对角矩阵，其中的对角元素是矩阵H的奇异值；

[0046] S53、分离出矩阵U中第1行到第L行，第1列到第P列的子矩阵为矩阵U1；同时分离出矩阵U中第2行到第L+1行，第1列到第P列的子矩阵为U2，以此得到U1和U2分别为：

[0047] U1＝U[1:L,1:P]

[0048] U2＝U[2:L+1,1:P]

[0049] 式中，P是分离后时频图谱确定的窄带干扰个数的2倍；

[0050] S54、构建矩阵Z＝[U1 U2]，然后对Z开展奇异值分解，得到

[0051]

[0052]

[0053] 式中：和 均是矩阵Z奇异值分解后的单位特征矩阵； 是矩阵Z奇异值分解后的对角矩阵；V11、V12、V21和V22是 的子矩阵；

[0054] S55、利用V12和V22构建矩阵ψ为

[0055] ψ＝‑V12V22‑1

[0056] S56、对ψ开展特征值分解，得到对应的特征值为μk(k＝1,2,…,P)，从而估计得到各窄带干扰分量的频率 为：

[0057]

[0058] 式中，arg(*)是计算复数的相角；

[0059] S57、利用最小二乘法计算各窄带干扰分量的幅值 和相角 具体计算公式为：

[0060] Y＝μc

[0061] Y＝[y(0)，y(1)，…，y(Ny‑1)]T

[0062]

[0063]

[0064]

[0065] 式中，c是中间变量矩阵，ck是c中元素。

[0066] 进一步的，所述信号帧内窄带干扰相位 通过相角 开展移相处理得到，其计算公式为：

[0067]

[0068] 式中，ΔT是信号帧和噪声帧的时延。

[0069] 综上所述，本发明与现有技术相比具有以下有益效果：本发明能有效分离出染噪PD信号时频谱图中窄带干扰和PD信号，具备良好的窄带干扰抑制效果，能有效恢复局部放电信号的波形特征，进而能提升窄带干扰影响下PD信号时频分析的准确率。

附图说明

[0070] 此处所说明的附图用来提供对本发明实施例的进一步理解，构成本申请的一部分，并不构成对本发明实施例的限定。在附图中：

[0071] 图1为PD信号的仿真波形图；

[0072] 图2为仿真PD波形的时频谱图；

[0073] 图3为染噪PD信号时频谱图中频率能量和时间的关系图；

[0074] 图4为分离后窄带干扰和原始PD信号的时频谱图；

[0075] 图5为仿真PD波形的窄带干扰抑制结果图；

[0076] 图6为混合干扰下的染噪PD波形图；

[0077] 图7为混合干扰下染噪PD信号的时频谱图分离结果图；

[0078] 图8为混合干扰下染噪PD波形的降噪结果图；

[0079] 图9为实测的染噪PD信号图；

[0080] 图10为实测PD信号的时频谱图；

[0081] 图11为实测PD信号的时频谱图分离结果图；

[0082] 图12为实测PD信号的窄带干扰抑制结果图；

[0083] 图13为本发明一个具体实施例的流程图。

具体实施方式

[0084] 为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施例和附图，对本发明作进一步的详细说明，本发明的示意性实施方式及其说明仅用于解释本发明，并不作为对本发明的限定。

[0085] 实施例：

[0086] 如图13所示，基于时频谱图分离的局部放电窄带干扰抑制方法，包括：S1、采样获取染噪PD信号；S2、采用广义S变换对染噪PD信号的离散信号进行时频变换处理，得到二维时频复数矩阵，其中，二维时频复数矩阵中横轴为时间采样点，纵轴为频率采样点；S3、对二维时频复数矩阵求模得到关于染噪PD信号的时频谱图；S4、采用数学形态学获取染噪PD信号时频谱图中各行分量的平滑分量和非平滑分量，将分离出的平滑分量和非平滑分量作为行分量构成新的矩阵，分离出窄带干扰和原始PD信号的时频谱图，以确定窄带干扰的数目与PD脉冲的时间区域；S5、将采样得到的染噪PD信号根据时间轴划定为含有PD信号的信号帧和仅含噪声的噪声帧，采用总体最小二乘不变旋转矢量技术分析噪声帧内时域数据，以提取窄带干扰的时域特征量；S6、重构出窄带干扰的时域波形，从染噪PD信号的时域波形中去除窄带干扰的时域波形，得到降噪后PD信号。

[0087] 本实施例中，信号x(t)的广义S变换结果G(τ,f,λ)被定义为：

[0088]

[0089] 式中，t和τ为两组时间变量；f是频率；λ是调节因子；w(t‑τ,f,λ)是高斯窗函数，对应的表达式为：

[0090]

[0091] 由于式(2)中窗函数的宽度和分析频率成反比，所以广义S变换在低频区域能获得更高的频率分辨率，而在高频区域能获得更高的时间分辨率，利于PD信号这类非平稳信号的时频分析。同时可以通过设定合适的λ，对广义S变换的最终时频变换结果进行人为调节。如果将λ设定在(0,1)范围内，那么可以增加式(2)中窗函数的宽度，达到提升频域分辨率的目的；如果将λ设定在(1,+∞)范围内，那么可以减小式(2)中窗函数的宽度，达到提升时域分辨率的目的。

[0092] 在实际使用广义S变换时，令x(n)是信号x的离散信号，同时定义f＝n/(NT)、τ＝iT，其中T是采样周期，N是x(n)的数据总量，得到x(n)的广义S变换结果为：

[0093]

[0094] 式中，i、m、n是x(n)的广义S变换结果中三组变量，被定义为0，1，…，N‑1。

[0095] 在利用式(3)对x(n)开展时频变换处理后，可以得到一个二维的时频复数矩阵，该矩阵的横轴是时间采样点，该矩阵的纵轴是频率采样点。对该二维时频复数矩阵进行求模处理，可以得到关于染噪PD信号的时频谱图。

[0096] 本实施例采用的广义S变换方法是一种基于S变换的改进时频变换方法，该方法在S变换方法的窗函数中加入了调节因子λ，不仅继承了S变换对非平稳信号的优异时频分析能力，同时增加了时频变换结果的人为可控性。本实施例采用的广义S变换方法同时具备小波变换方法的多分辨率分析和短时傅里叶变换方法的单频率分量分析的优点，并且不需要人为选择小波函数或窗函数的参数，减小了人为因素干扰。

[0097] 数学形态学中存在两种基本的变换形式：腐蚀和膨胀。染噪PD信号时频谱图中各行分量Gt(nGt)的膨胀和腐蚀变换结果分别表示为：

[0098]

[0099]

[0100] 式中，g(ig)是结构元素；⊕是膨胀运算符；⊙是腐蚀运算符；DGt是Gt(nGt)的定义域；Dg是g(ig)的定义域，nGt是数组Gt中某个索引，ig是数组g中某个索引。

[0101] 通过式(4)和式(5)的级联组合，可以得到数学形态学中开运算和闭运算分别为：

[0102]

[0103] (Gt·g)(nGt)＝((Gt⊕g)⊙g)(nGt) (7)

[0104] 式中，是开运算；·是闭运算。

[0105] 利用开、闭运算的组合可以对目标信号进行削峰平谷，从而分离出目标信号的平滑分量和非平滑分量。将式(6)和式(7)的开、闭运算进行级联组合，分别得到形态开‑闭滤波器运行结果[Foc(Gt)](nGt)和形态闭‑开滤波器运行结果[Fco(Gt)](nGt)为:

[0106]

[0107]

[0108] 在利用数学形态学算法对目标信号开展分量分离时，为了减小统计偏倚现象的影响，本实施例采用形态开‑闭滤波器和形态闭‑开滤波器的混合运算，得到信号Gt(nGt)中的平滑分量[F1(Gt)](nGt)和非平滑分量[F2(Gt)](nGt)分别为：

[0109] [F1(Gt)](nGt)＝([Foc(Gt)](nGt)+[Fco(Gt)](nGt))/2 (10)

[0110] [F2(Gt)](nGt)＝Gt(nGt)‑[F1(Gt)](nGt) (11)。

[0111] 本实施例基于数学形态学进行信号特征提取，利用合适的结构元素去采集目标信号的特征信息，从而刻画出目标信号的结构特征。

[0112] 本实施例采用总体最小二乘不变旋转矢量技术(total least squares‑estimation of signal parameters via rotational invariance techniques,TLE‑ESPRIT)计算窄带干扰的特征参数，其具体计算流程如下：

[0113] (1)将噪声帧内时域数据信号y(ny)构成Hankle矩阵H为：

[0114]

[0115] 式中，Ny是信号y(ny)的数据长度，L是H的列数，L的优选值为Ny/3。

[0116] (2)将H开展奇异值分解，得到

[0117] H＝UΣVT (13)

[0118] 式中，T是取矩阵的共轭转置；U和V分别是左、右正交矩阵；Σ是对角矩阵，其中的对角元素是矩阵H的奇异值。

[0119] (3)分离出矩阵U中第1行到第L行，第1列到第P列的子矩阵为矩阵U1；同时分离出矩阵U中第2行到第L+1行，第1列到第P列的子矩阵为U2，以此得到U1和U2分别为：

[0120] U1＝U[1:L,1:P] (14)

[0121] U2＝U[2:L+1,1:P] (15)

[0122] 式中，P是分离后时频图谱确定的窄带干扰个数的2倍。

[0123] (4)构建矩阵Z＝[U1 U2]，然后对Z开展奇异值分解，得到

[0124]

[0125]

[0126] 式中：和 均是矩阵Z奇异值分解后的单位特征矩阵； 是矩阵Z奇异值分解后的对角矩阵；V11、V12、V21和V22是 的子矩阵。

[0127] (5)利用V12和V22构建矩阵ψ为

[0128] ψ＝‑V12V22‑1 (18)

[0129] (6)对ψ开展特征值分解，得到对应的特征值为μk(k＝1,2,…,P)，从而估计得到各窄带干扰分量的频率 为：

[0130]

[0131] 式中，arg(*)是计算复数的相角。

[0132] (7)利用最小二乘法计算各窄带干扰分量的幅值 和相角 具体计算公式为：

[0133] Y＝μc (20)

[0134] Y＝[y(0)，y(1)，…，y(Ny‑1)]T (21)

[0135]

[0136]

[0137]

[0138] 式中，c是中间变量矩阵，ck是c中元素。

[0139] PD脉冲信号通常会呈现衰减振荡的波形特征，因此采用式(25)所示的单指数衰减振荡函数s1(t)和式(26)所示的双指数衰减振荡函数s2(t)对真实的局部放电信号进行模拟。

[0140]

[0141]

[0142] 式中，a1和a2是仿真PD信号的幅值；τ1和τ2是仿真PD信号的衰减常数；fs1和fs2是仿真PD信号的振荡频率。

[0143] 周期性窄带干扰snoise(t)的表达式为：

[0144]

[0145] 式中，Ak是窄带干扰的幅值；fk是窄带干扰的频率； 是窄带干扰的初始相角；kT是窄带干扰的数目。

[0146] 利用表1中参数构造4组仿真PD脉冲，如图1(a)所示。由于现场窄带干扰频带具有一定宽度，因此本实施例设计了多组不同频带的窄带干扰进行模拟测试。利用表2中参数构造5组窄带干扰，并将窄带干扰叠加在图1(a)中原始的PD仿真波形上，得到染噪的PD仿真波形如图1(b)所示，图1中波形的采样时长和采样频率分别设置为120μs和40MHz。值得说明的是，对比表1和表2可知，PD信号脉冲在脉冲序号1、3和窄带干扰的干扰序号3的频带存在重合现象，以模拟实际情况中局部放电信号频带与窄带干扰频带重合的情况。从图1中可以看出，在叠加了窄带干扰之后，原始PD信号几乎完全被窄带干扰淹没，因此无法直接从时域中将原始PD信号和窄带干扰进行分离，导致PD信号的识别和分析受到限制。

[0147] 表1仿真PD信号的参数表

[0148]脉冲序号 幅值/mV 衰减常数/μs 振荡频率/MHz 脉冲类型
1 7.5 1 2.5 s2(t)
2 1.5 1 4.2 s1(t)
3 1.5 0.8 2.5 s1(t)
4 7.5 0.8 4.2 s2(t)

[0149] 表2仿真窄带干扰的参数表

[0150] 干扰序号 幅值/mV 频率/MHz 初始相位/rad1 3 0.6 π/2
2 0.5 1.51 π/3
3 1 2.53 π/6
4 3 5.31 π/5
5 2 8.1 π/4

[0151] PD脉冲信号通常会呈现衰减振荡的波形特征，因此可以采用式(25)所示的单指数衰减振荡函数s1(t)和式(26)所示的双指数衰减振荡函数s2(t)对真实的局部放电信号进行模拟。

[0152] 利用广义S变换方法分别对图1(a)中原始的PD仿真波形和图1(b)中染噪的PD仿真波形进行时频分析，得到两者对应的时频谱图如图2所示。需要说明的是，本实施例将式(3)中λ优选设置为0.3。

[0153] 在图2的时频谱图中，原始PD信号和窄带干扰呈现出明显不同的特点。原始PD信号在时间轴方向上的跨度较小，同时该信号在频率轴方向上的跨度较大；窄带干扰在频率轴方向上的跨度较小，同时该信号在时间轴方向上的跨度较大。综上所述，可以利用原始PD信号和窄带干扰在时频谱图中的不同特征，对两者进行分离。进一步对图2(b)中染噪PD信号的时频谱图开展分析，可以发现：当窄带干扰和原始PD信号的时频能量存在混叠时，由于窄带干扰的能量较强，所以窄带干扰会掩盖原始PD信号的时频特征，难以直接在时频谱图中区分窄带干扰和PD信号。若直接在时频谱图中区分窄带干扰和PD信号，具有很大的局限性，难以准确划分出窄带干扰和PD信号对应的时频子区域。

[0154] 对图2(b)中染噪PD信号时频谱图的每行数据开展单独分析，即单独分析各频率能量和时间的变化关系曲线，得到2组典型的变化关系曲线如图3所示。其中图3(a)的频率能量仅由原始PD信号构成，图3(b)的频率能量由原始PD信号和窄带干扰共同构成。这里需要说明的是，由于染噪PD信号的首端和末端存在数据的截断效应，所以染噪PD信号的广义S变换模矩阵在时间序列的首端和末端都存在畸变，为了准确分析染噪PD信号广义S变换模矩阵的特征，本实施例在广义S变换模矩阵时间序列的首端和末端都删除一小部分采样点。

[0155] 分析图3可以看出，在各频率能量和时间的变化关系曲线中，窄带干扰会呈现为平滑分量，而原始PD信号是一个“丘状”结构，呈现为非平滑分量，本实施例利用该特点将两者进行区分，以分离出染噪PD信号时频谱图中窄带干扰与原始PD信号。本实施例采用数学形态学方法分离信号中平滑分量和非平滑分量，进而分离染噪PD信号时频谱图中窄带干扰和PD信号。

[0156] 通常局部放电的持续时间小于1μs，在测试回路中电容和电感的影响下，实际测试的PD波形的持续时间通常会大于1μs。为了有效分离染噪PD信号时频谱图中原始PD信号和窄带干扰，本实施例中结构元素选择为时间长度为5μs的扁平结构元素。

[0157] 依照本实施例，首先将染噪PD信号时频谱图中各行分量作为目标信号，利用数学形态学方法分离目标信号的平滑分量和非平滑分量，然后将分离出的平滑分量和非平滑分量作为行分量构成新的矩阵，以此分离出窄带干扰和原始PD信号的时频谱图，从而提取窄带干扰和原始PD信号的特征量。为了方便观察，本实施例对分离出的原始PD信号时频谱图进行了取模处理。

[0158] 利用数学形态学方法对图2(b)中染噪PD信号的时频谱图进行处理，分离出窄带干扰和原始PD信号的时频谱图分别如图4(a)和图4(b)所示。首先分析图4(a)中分离出的窄带干扰时频谱图，依照窄带干扰的时频特性，可以确定仿真的染噪PD信号中存在5组窄带干扰。接着分析图4(b)中分离出的原始PD信号时频谱图，由于部分PD信号和窄带干扰在染噪PD信号时频谱图中发生了时频能量的混叠，所以图4(b)中部分PD信号(仿真的局放脉冲1和局放脉冲3)的时频能量分布与图2(a)中原始PD信号的时频能量分布产生了差异，但是仍然具备原始PD信号的时频特性，因此依旧可以按照PD信号的时频特性，确定PD脉冲的时间区域。

[0159] 当采样的时间窗口较长时，PD信号仅会占据较小的时间，更多的时间是由窄带干扰等噪声所占据，所以可以将采样得到的染噪PD信号根据时间轴划定为含有PD信号的信号帧和仅含窄带干扰等噪声的噪声帧。

[0160] 由于染噪PD信号在噪声帧内不含有PD信号，因此可以利用TLE‑ESPRIT方法分析噪声帧内时域数据，以准确提取窄带干扰的时域特征量。需要说明的是，利用噪声帧内时域数据计算得到的相位仅对应噪声帧内的窄带干扰，为计算信号帧中相位以重构窄带干扰，需要将式(23)中 开展移相处理，得到信号帧内窄带干扰相位 为：

[0161]

[0162] 式中，ΔT是信号帧和噪声帧的时延。

[0163] 通过TLE‑ESPRIT提取窄带干扰的时域特征量后，然后利用式(28)可以重构出窄带干扰的时域波形，最后从染噪PD信号的时域波形中去除窄带干扰的时域波形，得到降噪后PD信号。

[0164] 对图1(b)所示含有窄带干扰的染噪PD信号利用本实施例进行降噪，得到降噪结果如图5(a)所示。为了便于比较，本实施例引入FFT谱最小熵解卷积滤波方法和自适应奇异值分解方法对仿真染噪PD信号进行窄带干扰降噪，分别得到对应的降噪结果如图5(b)和图5(c)所示。

[0165] 由图5可知，对于FFT谱最小熵解卷积滤波方法而言，由于傅里叶变换算法存在频谱泄露的特点，导致该方法的最终降噪结果存在明显的“边缘效应”(边缘处的噪声较大)，同时该方法中阈值法无法提取小幅值窄带干扰的频点，因此该方法难以抑制小幅值的窄带干扰，导致整个降噪波形中残留较多的窄带干扰能量。对于自适应奇异值分解方法而言，由于PD波形和小幅值窄带干扰的奇异值数值非常接近，甚至会导致奇异值的能量产生相互的混叠，所以最终的降噪结果存在窄带干扰抑制不干净的情况。利用本实施例对染噪的仿真PD波形进行降噪时，既可以最大程度地削弱窄带干扰的能量，同时也能保留PD信号的波形特征。

[0166] 为了对图5中各方法的实际降噪结果开展量化分析，本实施例采用了信噪比(signal to noise ratio，SNR)、波形相似系数(normalized correlation coefficient，NCC)和均方误差(mean square error，MSE)三组特征参数对仿真PD波形的降噪效果进行评估，计算得到图5中各方法降噪结果的评价参数如表3所示。对比表3中的各项参数可以看出，本实施例的SNR和NNC值最大，同时MSE值最小，说明本实施例对窄带干扰的降噪效果最好，对原始PD信号的还原效果最佳。

[0167] 表3降噪结果的评价参数表

[0168]

[0169] 通常情况下，实际采集的PD信号会同时含有白噪声干扰和窄带干扰，为了验证本实施例的适用性，先对图1(a)的仿真原始PD波形添加2dB的白噪声干扰，再加入仿真的窄带干扰，得到对应混合干扰下的PD波形如图6所示，此时计算得到染噪PD信号的信噪比为‑26.5915dB。

[0170] 利用本实施例分离图6中染噪PD波形的时频谱图，得到分离后窄带干扰和PD信号的时频谱图分别如图7(a)和图7(b)所示。从图7中可以看出，在染噪PD信号存在白噪声干扰的情况下，本实施例仍旧可以有效分离出窄带干扰和PD信号的时频谱图，其中分离出的窄带干扰时频谱图无异常变化，可以准确判断窄带干扰数量；分离出的PD信号时频谱图中会存在白噪声能量，但是PD信号的时频特征仍然清晰可见，可以准确划定PD信号的时频区域，确定染噪PD信号的信号帧和噪声帧。

[0171] 利用本实施例对图6中染噪PD信号的窄带干扰进行抑制，得到对应的降噪结果如图8所示。从图8中可以看出，此时已经可以明显观察到PD信号的波形，计算得到此时的信噪比已经降为1.4571dB。上述结果证明了：在存在白噪声干扰的情况下，本实施例仍旧可以有效抑制染噪PD信号中窄带干扰，从而提高染噪PD信号的信噪比。

[0172] 为了验证本实施例抑制实测染噪PD信号中窄带干扰的有效性，在实验室中搭建10kV交联聚乙烯电力电缆的工频局放测试系统，局放缺陷设置为电缆终端处存在半导电层突起缺陷，采样频率设置为100MHz，采样点数设置为4000。由于实验室近似于无噪环境，人为向实测PD信号中叠加5组窄带干扰，幅值分别设置为20、8、25、35、25mV，频率分别设置为
2.51、5.36、10.32、15.33、22.42MHz，相位分别设置为45°、60°、30°、30°和90°，得到实测的染噪PD信号如图9所示。从图9中可以看出，在窄带干扰的影响下，原始PD信号已经严重被污染，难以进行特征分析。

[0173] 利用广义S变换对图9中实测的染噪PD信号进行时频分析，得到对应的时频谱图如图10所示，从图10中可以看出，由于窄带干扰的能量较强，所以在实测染噪PD信号的时频谱图中，PD信号的时频特征已经完全被掩盖，直接在该时频谱图中划定窄带干扰和PD信号的方法难以实施，具备较大的局限性。

[0174] 利用本实施例对图10实测PD信号的时频谱图进行分离，得到分离后窄带干扰和PD信号对应的时频谱图如图11所示，从图11中可以看出，此时窄带干扰和PD信号的时频谱图得到了有效的分离，利用分离出的窄带干扰时频谱图可以确定窄带干扰个数为5，同时利用分离出的PD信号时频谱图可以划定信号帧和噪声帧，以确保利用TLS‑ESPRIT准确提取窄带干扰特征参数，重构窄带干扰，得到最终的降噪结果如图12(a)所示。

[0175] 为了说明本实施例的优越性，再分别利用FFT谱最小熵解卷积滤波法和自适应奇异值分解法对图9中实测的染噪PD信号进行降噪处理，分别得到降噪结果如图12(b)和图12(c)所示。从图12中可以看出，利用FFT谱最小熵解卷积滤波法进行降噪时，降噪结果中残余的窄带干扰能量较强，尤其是在信号的边缘部分；利用自适应奇异值分解法进行降噪时，虽然大部分窄带干扰能够被有效抑制，但是会残留部分的小幅值窄带干扰；而采用本实施例进行降噪时，降噪结果中窄带干扰的残余能量较小，同时PD信号的波形特征得到了保留。

[0176] 因为难以采集到完全没有噪声的PD信号，所以不能采用SNR、NCC、MSE等参数对实测PD信号的降噪效果进行量化分析，本实施例引入噪声抑制比ρ抑对图12中降噪效果进行评估。该参数可以表现降噪结果中PD信号的凸显程度，该值越大，说明降噪结果更好。ρ抑的具体定义为：

[0177]

[0178] 式中，σ1是降噪前信号的标准差，σ2是降噪后信号的标准差。

[0179] 计算得到本实施例的方法、FFT谱最小熵解卷积滤波法和自适应奇异值分解法的噪声抑制比结果分别为30.7559、7.3017、16.4680，通过对比可以发现，本实施例的ρ抑值更大，说明本实施例可以有效抑制窄带干扰，以恢复PD信号的波形特征。值得说明的是，由图11可知，在本实施例实测的染噪PD信号中，PD信号和窄带干扰存在明显的频带重叠，而本实施例依然可以取得较好的窄带干扰抑制效果。

[0180] 本实施例提出了一种基于时频谱图分离的局部放电窄带干扰抑制方法，该方法可以有效降低染噪PD信号中窄带干扰的能量，恢复PD信号的波形特征。本实施例的具体结论如下：

[0181] (1)通过分析染噪PD信号的时频谱图，明确当窄带干扰的能量较强时，窄带干扰会掩盖原始PD信号的时频特征，难以直接在染噪PD信号的时频谱图中辨识原始PD信号和窄带干扰。

[0182] (2)结合广义S变换和数学形态学，本实施例所提方法可以从染噪PD信号中分离出原始PD信号和窄带干扰信号的时频谱图，从而在时频域内准确辨识出窄带干扰和原始PD信号。

[0183] (3)在分离出原始PD信号和窄带干扰信号的时频谱图的基础上，利用TLS‑ESPRIT算法可以准确估计窄带干扰的特征参数，实现窄带干扰的抑制。

[0184] (4)仿真和实测染噪PD信号的降噪结果表明：对比传统的FFT谱最小熵解卷积滤波法和自适应奇异值分解法，本实施例具备更好的窄带干扰抑制效果，有利于小幅值窄带干扰的抑制。

[0185] 利用本实施例的方法、FFT谱最小熵解卷积滤波方法和自适应奇异值分解方法对仿真和实测的染噪PD信号开展降噪处理，通过对比降噪结果，验证了本实施例可以有效抑制染噪PD信号中窄带干扰，并且降噪结果中残余的窄带干扰能量较小。

[0186] 以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。