基于廣義S變換和隨機(jī)子空間的局放窄帶干擾抑制方法

宋立業(yè)， 蒲霄祥， 李希桐

(遼寧工程技術(shù)大學(xué)電氣與控制工程學(xué)院， 遼寧 葫蘆島 125105)

1 引言

局部放電(Partial Discharge，PD)監(jiān)測(cè)是有效診斷電力電纜絕緣狀態(tài)的手段之一，一般情況下PD的信號(hào)能量較弱，會(huì)被現(xiàn)場(chǎng)的電磁干擾所掩蓋，阻礙了PD信號(hào)的提取識(shí)別[1,2]。PD信號(hào)面臨的電磁干擾主要分為隨機(jī)脈沖干擾、周期性窄帶干擾和白噪聲三大類，其中周期性窄帶干擾擁有噪聲能量強(qiáng)和持續(xù)時(shí)間長(zhǎng)的特點(diǎn)，會(huì)嚴(yán)重污染PD信號(hào)，因此研究PD信號(hào)中窄帶干擾抑制方法具有重大的意義[3]。

目前，國(guó)內(nèi)外學(xué)者針對(duì)PD信號(hào)的窄帶干擾抑制方法開展了大量研究。文獻(xiàn)[4]提出利用快速傅里葉變換配合頻域中閾值抑制窄帶干擾，該方法操作簡(jiǎn)單、計(jì)算速度快，但是容易受到頻譜泄露的影響造成嚴(yán)重的邊緣效應(yīng)。文獻(xiàn)[5,6]提出利用小波分解方法開展PD信號(hào)去噪，該方法可以有效分析PD信號(hào)的時(shí)頻特征，進(jìn)而取得較好的去噪效果，但是該方法需要人為設(shè)定小波基、分解層數(shù)和各層閾值，去噪效果受人為因素影響較大。文獻(xiàn)[7]提出利用經(jīng)驗(yàn)?zāi)B(tài)分解方法對(duì)PD信號(hào)進(jìn)行自適應(yīng)分解去噪，該方法不需要選擇基函數(shù)，能更好地適應(yīng)多種PD信號(hào)，但是該方法存在端點(diǎn)效應(yīng)和模態(tài)混疊的問題。文獻(xiàn)[8,9]提出利用奇異值分解實(shí)現(xiàn)窄帶干擾和PD信號(hào)的分離，該方法僅需要確定奇異值閾值便可以重構(gòu)出純凈的PD信號(hào)，但是該方法難以準(zhǔn)確給出奇異值閾值，同時(shí)該方法在窄帶干擾幅值較低以及干擾和PD信號(hào)的頻率出現(xiàn)混疊時(shí)去噪效果較差。

廣義S變換具有良好的時(shí)頻分辨能力，能有效地分離出時(shí)頻特征不同的PD信號(hào)和窄帶干擾，因此被逐漸用于PD信號(hào)的窄帶干擾抑制中。文獻(xiàn)[10]提出利用廣義S變換模時(shí)頻矩陣的能量分布重構(gòu)窄帶干擾進(jìn)行去噪，該方法可以有效在時(shí)頻域中分離出PD信號(hào)和窄帶干擾，但是沒有考慮窄帶干擾相位的影響。文獻(xiàn)[11]提出利用廣義S變換配合奇異值分解算法實(shí)現(xiàn)窄帶干擾抑制，雖然解決了奇異值分解算法中奇異值閾值的選取問題，但是仍在窄帶干擾和PD信號(hào)的頻率出現(xiàn)混疊時(shí)去噪效果較差。

綜上所述，本文提出一種基于廣義S變換和隨機(jī)子空間的局部放電窄帶干擾抑制方法。該方法首先通過(guò)廣義S變換得到染噪PD信號(hào)的模時(shí)頻矩陣，然后根據(jù)窄帶干擾和PD信號(hào)各自的時(shí)頻特征確定窄帶干擾數(shù)量和分離出不含PD的時(shí)間片段，接著利用隨機(jī)子空間算法進(jìn)行窄帶干擾數(shù)據(jù)重構(gòu)，最后提取出純凈的PD信號(hào)。分別利用本文方法、廣義S變換模矩陣方法和頻率切片小波變換方法對(duì)仿真和實(shí)測(cè)的染噪PD數(shù)據(jù)進(jìn)行窄帶干擾抑制，對(duì)比結(jié)果顯示，相比于傳統(tǒng)方法，本文方法能更好地抑制局部放電信號(hào)中的窄帶干擾。

2 算法原理

2.1 廣義S變換

S時(shí)頻變換早期由Stockwell提出，該方法集合了短時(shí)傅里葉變換和連續(xù)小波變換的優(yōu)點(diǎn)[10,11]。S變換中窗函數(shù)選用了參數(shù)可變的高斯函數(shù)，該函數(shù)的時(shí)寬和頻率成反比，幅值與頻率成正比，因此S變換可以在低頻區(qū)間獲得較好的頻率分辨率，而在高頻區(qū)間獲得較好的時(shí)間分辨率，擁有優(yōu)異的時(shí)頻分析能力，對(duì)于提取PD信號(hào)這類非平穩(wěn)信號(hào)的時(shí)頻分布特征具有良好的效果。

信號(hào)y(t)的S變換可以表示為:

(1)

式中,f為頻率；ω(t-τ,f)為高斯窗函數(shù)；t、τ為時(shí)間變量。

S變換中窗函數(shù)表達(dá)式為：

(2)

從式(2)中可以看出，S變換中窗函數(shù)的幅值和時(shí)寬是隨頻率變化而變化的，因此S變換中時(shí)頻分辨率能夠跟隨頻率變化而變化，能改善短時(shí)傅里葉變換中不變的時(shí)頻分辨率，同時(shí)不需要考慮連續(xù)小波變換中基函數(shù)的選擇問題，憑借上述優(yōu)點(diǎn)，S變換具有良好的使用前景。

為了使S變換能夠適用更多的使用場(chǎng)景，Pinnegar設(shè)計(jì)了調(diào)節(jié)因子λ(λ>0)對(duì)S變換中高斯窗函數(shù)進(jìn)行了改進(jìn)[10,11]，得到新的窗函數(shù)為：

(3)

使用式(3)作為窗函數(shù)的S變換被稱為廣義S變換，從式(3)中可以看出，廣義S變換可以通過(guò)調(diào)節(jié)λ實(shí)現(xiàn)高斯窗函數(shù)和f的變化率調(diào)節(jié)。當(dāng)λ<1時(shí)，高斯窗的時(shí)寬增加，幅值減少，廣義S變換的頻率分辨率增加，時(shí)間分辨率下降；當(dāng)λ>1時(shí)，高斯窗的時(shí)寬減少，幅值增加，廣義S變換的頻率分辨率下降，時(shí)間分辨率增加，由此實(shí)現(xiàn)時(shí)頻分辨率的有效調(diào)節(jié)。

由于實(shí)際信號(hào)通常為離散數(shù)據(jù)，因此需要對(duì)廣義S變換做離散化處理，令f=n/(NT)和τ=iT，其中T為采樣周期，N為采樣個(gè)數(shù)，得到離散的廣義S變換為：

(4)

式中，n，i=0，1，…，N-1；Y是y的傅里葉系數(shù)。

通過(guò)式(4)得到信號(hào)的廣義S變換復(fù)數(shù)矩陣。為了能直觀分析信號(hào)在時(shí)頻域中的能量分布，對(duì)該復(fù)數(shù)矩陣進(jìn)行求模，得到廣義S變換模矩陣(Generalized S-transform Modular Matrix，GSMM)。該矩陣能直接反映信號(hào)能量的時(shí)間和頻率分布特性，因此可以利用GSMM對(duì)信號(hào)的時(shí)頻特征進(jìn)行深度分析。

2.2 隨機(jī)子空間

隨機(jī)子空間算法可以有效地計(jì)算離散系統(tǒng)的模態(tài)參數(shù)，并且有著計(jì)算精度高、受采樣條件影響小和步驟簡(jiǎn)單的優(yōu)點(diǎn)[12]。其中基于協(xié)方差驅(qū)動(dòng)的隨機(jī)子空間算法憑借著優(yōu)異的計(jì)算效率被廣泛使用[13]，本文將其引入用于窄帶干擾參數(shù)估計(jì)。

將離散信號(hào)yk以離散系統(tǒng)的狀態(tài)方程可以表示為：

(5)

式中，A為離散系統(tǒng)的狀態(tài)矩陣；C為離散系統(tǒng)的輸入矩陣；xk和yk分別為離散系統(tǒng)k時(shí)刻的狀態(tài)量和輸出量；wk為離散系統(tǒng)中白噪聲；vk為測(cè)量中白噪聲。對(duì)應(yīng)本文中，yk是離散的窄帶干擾信號(hào)。

將yk建立三個(gè)時(shí)間矩陣分別為：

(6)

(7)

(8)

式中，Yp為“過(guò)去”時(shí)間矩陣；Yf1為第一個(gè)“未來(lái)”時(shí)間矩陣；Yf2為第二個(gè)“未來(lái)”時(shí)間矩陣；b為時(shí)間矩陣的列數(shù)，該值越大，說(shuō)明離散系統(tǒng)輸出數(shù)據(jù)量越多，此時(shí)各參數(shù)估計(jì)越精準(zhǔn)，但是yk的數(shù)據(jù)量是有限的，因此本文取b>10(N-b-1)，得到b為：

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

(9)

式中，F(xiàn)loor是向下取整。

將三個(gè)時(shí)間矩陣Yp、Yf1和Yf2構(gòu)建托普利茨矩陣T1和T2為：

T1=Yf1YpT

(10)

T2=Yf2YpT

(11)

接著將T1開展奇異值分解為：

(12)

式中，S1、S2分別為信號(hào)和噪聲主導(dǎo)的奇異值對(duì)角矩陣；U1、U2分別為信號(hào)和噪聲主導(dǎo)的左單位奇異矩陣；V1、V2分別為信號(hào)和噪聲主導(dǎo)的右單位奇異矩陣。

從統(tǒng)計(jì)理論上分析，wk和vk應(yīng)該是無(wú)相關(guān)性的兩組噪聲，由此根據(jù)離散系統(tǒng)的理論可得：

(13)

式中，Oa為觀測(cè)矩陣；G=E(x(k+1)ykT)，E為期望；Гa為控制矩陣；a=N-b-1。

將式(12)和式(13)進(jìn)行聯(lián)合分析可得：

(14)

再依據(jù)式(13)可得：

T2=OaAΓa

(15)

進(jìn)一步聯(lián)合式(14)和式(15)得到A為：

(16)

對(duì)A開展特征值分解為：

A=ψΛψ-1

(17)

式中，Λ=diag(zs)，s=1，2，…，p；zs是A的第s個(gè)特征值；p為式(5)中系統(tǒng)的階數(shù)，對(duì)應(yīng)本文為窄帶干擾個(gè)數(shù)的2倍；ψ為A的特征矩陣。

利用最小二乘法求解下式即可得到參數(shù)Bs。

(18)

利用zs和Bs得到窄帶干擾的幅值Qs、頻率fs和相位θs分別為：

Qs=2|Bs|

(19)

(20)

(21)

由于窄帶干擾為實(shí)信號(hào)，因此確定fs>0的對(duì)應(yīng)數(shù)據(jù)組為真實(shí)數(shù)據(jù)，進(jìn)而估計(jì)得到窄帶干擾的幅值、頻率和相位。

2.3 PD窄帶抑制流程

本文結(jié)合廣義S變換和隨機(jī)子空間提出PD窄帶干擾抑制方法的具體步驟為：

(1)將染噪PD信號(hào)開展廣義S變換，得到對(duì)應(yīng)的廣義S變換模矩陣GSMM。

(2)借助窄帶干擾時(shí)間分布長(zhǎng)，頻率能量集中的特點(diǎn)確定窄帶干擾個(gè)數(shù)。

(3)借助PD信號(hào)時(shí)間分布短，頻率能量分布寬的特點(diǎn)確定染噪PD信號(hào)中無(wú)PD的最長(zhǎng)時(shí)間片段。

(4)利用隨機(jī)子空間算法對(duì)染噪PD信號(hào)中無(wú)PD的最長(zhǎng)時(shí)間片段進(jìn)行處理，估計(jì)窄帶干擾參數(shù)。

(5)利用估計(jì)的窄帶干擾參數(shù)值對(duì)窄帶干擾信號(hào)進(jìn)行重構(gòu)，在時(shí)域中將染噪PD信號(hào)減去窄帶干擾的重構(gòu)值，得到窄帶干擾抑制后的PD信號(hào)。

3 仿真測(cè)試

3.1 仿真局放模型

PD信號(hào)通常呈現(xiàn)衰減振蕩的特性，因此可以用如式(22)和式(23)所示的單指數(shù)衰減振蕩型和雙指數(shù)衰減振蕩型信號(hào)模型進(jìn)行模擬[14]。

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

(22)

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

(23)

式中，C為脈沖幅值；η為衰減系數(shù)；fg為振蕩頻率。

本節(jié)模擬4個(gè)PD脈沖信號(hào)，各脈沖信號(hào)參數(shù)如表1所示，其中脈沖Ⅰ和Ⅲ是單指數(shù)衰減振蕩型脈沖，脈沖Ⅱ和Ⅳ是雙指數(shù)衰減振蕩型脈沖，采樣頻率為30 MHz。

表1 仿真PD脈沖參數(shù)

仿真中窄帶干擾由4個(gè)不同頻率、幅值和初始相位的正弦波疊加組成，根據(jù)文獻(xiàn)[11、15]的仿真窄帶干擾參數(shù)范圍，本文將仿真中窄帶干擾的參數(shù)值設(shè)置如表2所示。由此得到純凈的PD信號(hào)、窄帶干擾信號(hào)和染噪PD信號(hào)如圖1所示。從圖1中可以看出，在染噪PD信號(hào)中，由于窄帶干擾的存在，PD信號(hào)幾乎完全被淹沒，難以直接進(jìn)行識(shí)別和提取。

表2 仿真窄帶干擾參數(shù)

圖1 仿真的PD信號(hào)

3.2 局放去噪

由于廣義S變換仍屬于短時(shí)傅里葉變換，因此其時(shí)間分辨率和頻率分辨率并不能同時(shí)達(dá)到最高，本文將λ取為0.4以同時(shí)獲得較好的時(shí)間分辨率和頻率分辨率。通過(guò)廣義S變換處理圖1(c)中的染噪PD信號(hào)，得到染噪PD信號(hào)的GSMM如圖2所示。從圖2中可以看出窄帶干擾信號(hào)和PD信號(hào)存在明顯不同的特征：窄帶干擾信號(hào)時(shí)間上分布較長(zhǎng)，頻率上分布較集中；PD信號(hào)時(shí)間上分布較短，頻率上分布較寬。

以此可以通過(guò)分析圖2中染噪PD信號(hào)的GSMM確定窄帶干擾的個(gè)數(shù)為4，同時(shí)可以在時(shí)間橫軸上進(jìn)行區(qū)域劃分，得到無(wú)PD區(qū)域。對(duì)于圖2而言，橫軸上采樣點(diǎn)數(shù)為1 400～2 400區(qū)域可以被視為無(wú)PD區(qū)域，該區(qū)域中主導(dǎo)信號(hào)為窄帶干擾信號(hào)。將最長(zhǎng)時(shí)段的無(wú)PD區(qū)域作為窄帶干擾參數(shù)估計(jì)時(shí)間段，選擇圖1(c)中染噪PD信號(hào)中該時(shí)間段的數(shù)據(jù)進(jìn)行隨機(jī)子空間算法處理，得到窄帶干擾參數(shù)估計(jì)值如表3所示。對(duì)比表2和表3可以看出，該方法可以有效地估計(jì)窄帶干擾的各參數(shù)值。值得說(shuō)明的是，由于僅截取了染噪PD信號(hào)中部分時(shí)段區(qū)段進(jìn)行參數(shù)估計(jì)，因此估計(jì)得到的相位和原始參數(shù)的相位是不一致的。

表3 仿真染噪PD信號(hào)的窄帶干擾參數(shù)估計(jì)值

利用表3中窄帶干擾參數(shù)估計(jì)值重構(gòu)整個(gè)時(shí)間段的窄帶干擾，重構(gòu)的窄帶干擾Dg如下：

(24)

式中，t1為截取區(qū)段的起始時(shí)刻。

將染噪PD信號(hào)減去重構(gòu)的窄帶干擾得到去噪結(jié)果如圖3(a)所示，為了進(jìn)行對(duì)比，利用廣義S變換模矩陣方法[10]和頻率切片小波變換方法[16]處理圖1(c)中染噪PD信號(hào)，得到去噪結(jié)果分別如圖3(b)和圖3(c)所示。從圖3的對(duì)比結(jié)果可以看出，對(duì)于廣義S變換模矩陣方法而言，該方法沒有考慮窄帶干擾相位的影響，因此對(duì)于本文中窄帶干擾相位不為0的情況，該方法無(wú)法實(shí)現(xiàn)窄帶干擾抑制；對(duì)于頻率切片小波變換方法而言，去噪后波形存在較大殘余噪聲，同時(shí)波形中存在明顯的邊緣效應(yīng)，整體去噪效果較差；本文方法去噪后波形的噪聲較小，PD信號(hào)的細(xì)節(jié)得到保留，利于后續(xù)PD信號(hào)的提取識(shí)別分析。

圖3 仿真染噪PD信號(hào)的窄帶干擾抑制結(jié)果

為了進(jìn)一步說(shuō)明3種方法的去噪效果，本文引入信噪比SNR、均方誤差MSE和波形相似度NCC三個(gè)評(píng)價(jià)參數(shù)開展分析[17]，其具體計(jì)算公式為：

(25)

(26)

(27)

式中，g(n)為原始PD信號(hào)；d(n)為各方法去噪后的信號(hào)。其中SNR越大，MSE越小，NNC越接近于1時(shí)，去噪效果越好。

計(jì)算得到圖3中各方法去噪后波形的SNR、MSE和NCC如表4所示，從表4中可以看出，和廣義S變換模矩陣方法、頻率切片小波變換方法相比，本文方法能更好地去除染噪PD信號(hào)中窄帶干擾，去噪后波形更好地保留了原始PD信號(hào)的波形特征，利于后續(xù)的分析研究。

表4 窄帶干擾抑制結(jié)果對(duì)比

4 實(shí)測(cè)研究

為了驗(yàn)證本文方法在實(shí)際測(cè)試中的可行性，在實(shí)驗(yàn)室中對(duì)10 kV電纜開展工頻局放測(cè)試，電纜終端頭中制作有半導(dǎo)電層突刺缺陷[18]，監(jiān)測(cè)方法為高頻電流法[7]，采樣率設(shè)置為200 MHz，其具體接線圖如圖4所示。

圖4 工頻局放測(cè)試接線圖

由于實(shí)驗(yàn)室中采集的PD信號(hào)噪聲較小，因此在采集的PD信號(hào)中添加多個(gè)窄帶干擾信號(hào)。文獻(xiàn)[19]指出高頻電流法受到的窄帶干擾主要包括中波段0.5～1.6 MHz、短波段 2.3～25 MHz和調(diào)頻段88～108 MHz的廣播信號(hào)。因此本文在此基礎(chǔ)上隨機(jī)選擇窄帶干擾信號(hào)的幅值分別為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，得到實(shí)測(cè)的染噪PD信號(hào)如圖5所示，從圖5中可以看出，由于窄帶干擾的污染，PD信號(hào)無(wú)法直接被識(shí)別提取。

圖5 實(shí)測(cè)的染噪PD信號(hào)

利用本文方法對(duì)圖5中染噪PD信號(hào)進(jìn)行去噪處理，得到染噪PD信號(hào)的GSMM圖如圖6所示，借助PD信號(hào)和窄帶干擾的時(shí)頻特征，可以確定圖6中窄帶干擾數(shù)目為4，同時(shí)很容易在時(shí)間橫軸上確定無(wú)PD區(qū)域。對(duì)于圖6而言，時(shí)間橫軸上采樣點(diǎn)數(shù)為600～1 400區(qū)域可以被視為無(wú)PD區(qū)域，將圖5中染噪PD信號(hào)對(duì)應(yīng)時(shí)間段的數(shù)據(jù)利用隨機(jī)子空間算法進(jìn)行窄帶干擾參數(shù)估計(jì)，得到窄帶干擾參數(shù)估計(jì)結(jié)果如表5所示，利用該結(jié)果對(duì)窄帶干擾進(jìn)行重構(gòu)，得到去噪后波形如圖7(a)所示。為了進(jìn)一步說(shuō)明本文方法的優(yōu)越性，再利用廣義S變換模矩陣方法和頻率切片小波變換方法對(duì)圖5中染噪PD信號(hào)進(jìn)行去噪，得到對(duì)應(yīng)去噪結(jié)果如圖7(b)和圖7(c)所示。從圖7的對(duì)比結(jié)果中可以看出，廣義S變換模矩陣方法抑制窄帶干擾失敗，沒有提取出局放波形；頻率切片小波方法的去噪結(jié)果中存在明顯的窄帶干擾抑制不干凈現(xiàn)象，同時(shí)存在明顯的邊緣效應(yīng)，因此上述2種傳統(tǒng)方法的去噪效果較差。相對(duì)于傳統(tǒng)方法而言，本文方法可以更加有效地提取出PD波形，殘余噪聲更小。

圖6 實(shí)測(cè)染噪PD信號(hào)的GSMM圖

表5 實(shí)測(cè)染噪PD信號(hào)的窄帶干擾參數(shù)估計(jì)值

圖7 實(shí)測(cè)染噪PD信號(hào)的窄帶干擾抑制結(jié)果

為了對(duì)圖7中降噪結(jié)果進(jìn)行量化分析，本文引入了噪聲抑制比ρ[9]如式(28)所示。噪聲抑制比可以顯示出窄帶干擾抑制前后有效信號(hào)的凸顯程度，該值越大，說(shuō)明窄帶干擾抑制結(jié)果越好。

(28)

式中，σ1和σ2分別為窄帶干擾抑制前、后的信號(hào)標(biāo)準(zhǔn)差。

通過(guò)計(jì)算得到圖7中3種方法的噪聲抑制比分別為10.541 1、-2.854 8和8.188 0，可見，本文方法的噪聲抑制比最大，對(duì)窄帶干擾的抑制結(jié)果最好。

5 結(jié)論

(1)廣義S變換能夠?qū)⑷驹隤D信號(hào)從時(shí)域轉(zhuǎn)化到時(shí)頻域中，并具有較好的時(shí)頻分辨能力，借助PD脈沖和窄帶干擾的不同時(shí)頻特征可以確定染噪PD信號(hào)中窄帶干擾數(shù)目和無(wú)PD時(shí)間片段。

(2)利用隨機(jī)子空間算法和窄帶干擾數(shù)目處理無(wú)PD時(shí)間片段可以精確估計(jì)窄帶干擾參數(shù)，進(jìn)而有效重構(gòu)染噪PD信號(hào)中窄帶干擾，對(duì)染噪PD信號(hào)實(shí)現(xiàn)窄帶干擾抑制。

(3)仿真和實(shí)測(cè)結(jié)果表明，相比于廣義S變換模矩陣方法和頻率切片小波變換方法，本文所提方法能夠有效抑制染噪PD信號(hào)中的窄帶干擾，去噪后殘余噪聲較小，PD波形恢復(fù)效果更好。