結(jié)合廣義S變換和快速獨立分量分析的局放信號中窄帶干擾抑制方法

李向群，伍亞萍，莊旭菲，吉飛敏，張燕

（1. 西北民族大學(xué) 電氣工程學(xué)院, 甘肅省 蘭州市 730030；2. 內(nèi)蒙古工業(yè)大學(xué) 信息工程學(xué)院,內(nèi)蒙古自治區(qū) 呼和浩特市 010051）

0 引言

電力電纜是輸配電系統(tǒng)的重要組成部分[1]，而局部放電（partial discharge，PD）是電力電纜絕緣劣化的早期表現(xiàn)形式之一，因此PD信號的監(jiān)測和分析已經(jīng)成為電力電纜絕緣評估的有效手段[2]。但是，由于實際采集的PD信號能量較小，而實際現(xiàn)場的電磁干擾能量較大，所以PD信號通常會被電磁干擾淹沒，嚴(yán)重妨礙PD信號的監(jiān)測和分析。常見的電磁干擾主要為隨機脈沖干擾[3]、白噪聲[4-5]和周期性窄帶干擾[6-7]。其中周期性窄帶干擾有著較強的能量和較長的持續(xù)時間，通常會將采集的PD信號完全淹沒，因此PD信號中周期性窄帶干擾的抑制技術(shù)具有重大的研究意義[8]。

在PD信號中周期性窄帶干擾的抑制技術(shù)上，國內(nèi)外學(xué)者已經(jīng)開展了大量的研究工作。文獻[8-9]利用快速傅里葉變換（fast Fourier transform，F(xiàn)FT）方法有效削弱了周期性窄帶干擾的能量，但是由于頻譜泄露問題的存在，所以該類方法難以獲得較好的干擾抑制效果。文獻[10]利用小波分解方法的優(yōu)異時頻分析能力提升了干擾抑制效果，但是由于PD信號波形存在多樣性的特點，所以該類方法難以選擇合適的小波基函數(shù)。文獻[11]利用經(jīng)驗?zāi)B(tài)分解方法的自適應(yīng)分解能力對窄帶干擾進行抑制，避免了基函數(shù)的選取問題，但是由于端點效應(yīng)和模態(tài)混疊問題的存在，所以該類方法的運用范圍極其有限。文獻[12-13]利用奇異值分解方法將信號分為窄帶干擾子空間和PD信號子空間，該類方法利用奇異值分解方法的信號自適應(yīng)分解能力實現(xiàn)了窄帶干擾的抑制，但是該類方法存在奇異值閾值選取困難和部分窄帶干擾抑制效果差的問題。

為了解決現(xiàn)有技術(shù)存在的問題，本文提出基于廣義S變換和快速獨立分量分析的局放信號窄帶干擾抑制方法。該方法首先利用廣義S變換獲取窄帶干擾的特征區(qū)域，然后在窄帶干擾的特征區(qū)域中，利用Candan算法對窄帶干擾頻率進行精準(zhǔn)估計，最后利用改進的快速獨立分量分析方法分離PD信號，從而抑制PD信號中窄帶干擾。仿真和實際的測試結(jié)果說明，相比于傳統(tǒng)方法，本文方法能夠有效提升窄帶干擾的抑制效果。

1 窄帶干擾抑制方法

1.1 仿真信號

大量的研究結(jié)果表明，PD信號可以采用單指數(shù)振蕩衰減函數(shù)和雙指數(shù)振蕩衰減函數(shù)進行模擬分析研究[14]，即仿真PD信號的函數(shù)表達式為

式中：fp是PD信號的振蕩頻率；τ是PD信號的衰減系數(shù)；A是PD信號的幅值系數(shù)。

通過式（1）和式（2）的函數(shù)模型和表1中仿真PD脈沖的參數(shù)值，得到仿真PD脈沖如圖1（a）所示，其中脈沖1和脈沖2為單指數(shù)振蕩衰減函數(shù)模型，脈沖3和脈沖4為雙指數(shù)振蕩衰減函數(shù)模型。該仿真中采樣時長設(shè)置為50 μs，采樣頻率fs設(shè)置為50 MHz。

表1 仿真PD脈沖的參數(shù)Table 1 Parameters of simulated PD pulse

圖1 仿真的PD信號Fig. 1 Simulated PD signal

一般情況下，周期性窄帶干擾會表現(xiàn)為正弦信號或余弦信號，因此通常可以采用三角周期函數(shù)進行模擬[15]，具體的函數(shù)表達式為

式中：m是窄帶干擾的總數(shù)；Bi是窄帶干擾的幅值；fni是窄帶干擾的頻率；θi是窄帶干擾的相位。

本文將m設(shè)置為3，并且各仿真窄帶干擾的設(shè)置值如表2所示，以此得到仿真的窄帶干擾波形如圖1（b）所示。并將該窄帶干擾疊加到仿真PD脈沖中，得到染噪PD信號波形如圖1（c）所示，從圖1（c）中可以看出，原始的PD信號已經(jīng)完全被窄帶干擾所掩蓋，無法進行識別分析，因此需要進行窄帶干擾抑制。

表2 仿真窄帶干擾的參數(shù)Table 2 Parameters of simulated narrowband interference

1.2 廣義S變換

S變換是由Stockwell研究得到的一種信號時頻分析方法，該方法同時兼顧了短時傅里葉變換和小波變換的優(yōu)點[16-17]。廣義S變換在S變換的基礎(chǔ)上引入了調(diào)節(jié)因子λ，以此實現(xiàn)時頻分辨率的人為調(diào)節(jié)，獲得更加適宜的時頻分辨率[18]。

定義廣義S變換的窗函數(shù)w(η-t,f,λ)為

從式（4）中可以看出，通過人為設(shè)定λ值可以改變窗口長度，當(dāng)λ>1時，可以減小窗口長度，從而提升時域分辨率；當(dāng)λ<1時，可以增加窗口長度，從而提升頻率分辨率；當(dāng)λ=1時，廣義S變換退化為S變換。

對于信號x的離散形式x(i)，i=0，1，…，N-1而言，得到對應(yīng)的廣義S變換為

式中：T是采樣周期；n=0，1，…，N-1；N是采樣總長度。

通過式（5）可以得到x(i)的廣義S變換時頻分析矩陣，該矩陣中行向量代表時間；列向量代表頻率；矩陣中元素為復(fù)數(shù)，包含了廣義S變換結(jié)果的幅值信息和相位信息。為了方便對廣義S變換結(jié)果進行分析，對廣義S變換時頻分析矩陣進行求模處理得到模矩陣為

通過對式（6）中模矩陣進行分析，可以獲取信號中各頻率分量的時間分布情況。

1.3 窄帶干擾特征量分析

通過廣義S變換可以對染噪PD信號進行時頻分析，從而在時頻分布圖譜上對窄帶干擾和PD信號進行單獨分析。由于S變換的時頻分辨率完全由分析頻率所決定，無法進行人為調(diào)節(jié)，所以將廣義S變換的λ設(shè)置為0.4，以獲得合適的時間分辨率和頻率分辨率，便于在時頻分布圖譜中分析窄帶干擾和PD信號的時頻特征，從而提取出窄帶干擾的特征區(qū)域。

圖2為圖1（c）中染噪PD信號的廣義S變換模矩陣，即染噪PD信號的時頻分布圖譜。通過對染噪PD信號的時頻分布圖譜進行分析，可以劃分出窄帶干擾的特征區(qū)域，從而準(zhǔn)確提取窄帶干擾的特征量，實現(xiàn)窄帶干擾抑制。在染噪PD信號的時頻分布圖譜中，周期性窄帶干擾的時間分布較長，頻率分布較集中；PD信號的時間分布較短，頻率分布較廣，因此可以很容易地在時頻分布圖譜中區(qū)分出窄帶干擾和PD信號。在區(qū)分出窄帶干擾和PD信號后，可以確定窄帶干擾的總數(shù)為3，然后對各窄帶干擾對應(yīng)的時頻矩陣區(qū)域進行分析，將沒有和PD信號發(fā)生混疊的時頻子矩陣區(qū)域劃出，作為各窄帶干擾的特征區(qū)域，以此得到窄帶干擾特征區(qū)域Ⅰ、Ⅱ和Ⅲ。

圖2 仿真PD信號的時頻分布Fig. 2 Time frequency distribution picture of simulated PD signal

為了在各窄帶干擾特征區(qū)域中準(zhǔn)確估計窄帶干擾的頻率，本文引入了Candan算法[19]。該算法在傅里葉頻率估計算法的基礎(chǔ)上，考慮了原有算法存在的柵欄效應(yīng)問題，能取得較高的頻率估計精度，同時該算法避免了加窗賦權(quán)和多次迭代的計算步驟，具有較低的計算復(fù)雜度，利于工程中實現(xiàn)，其具體計算步驟如下：

1）對窄帶干擾特征區(qū)域進行邊緣標(biāo)記，得到窄帶干擾特征區(qū)域的頻率范圍為[fmin，fmax]，采樣點范圍為[nmin，nmax]，得到采樣點范圍內(nèi)長度N1=nmax-nmin+1，同時獲得離散頻率點范圍為[Floor(N1·fmin/fs)，Ceil(N1·fmax/fs)]，F(xiàn)loor(*)是指向下取整，Ceil(*)是指向上取整，以此得到窄帶干擾的特征采樣點內(nèi)數(shù)據(jù)為y(i)，i=0，1，…，N1-1，并且y(i)=x(i+nmin)。

2）對y(i)中離散頻率點范圍[Floor(N1·fmin/fs)，Ceil(N1·fmax/fs)]進行離散傅里葉變換，得到傅里葉系數(shù)為

3）確定Y(k)的幅值最大值位置，并記錄對應(yīng)的離散頻率點為p，對應(yīng)的傅里葉系數(shù)為Y(p)，同時確定附近的傅里葉系數(shù)分別為Y(p-1)和Y(p+1)。

4）計算頻率修正值為

5）確定窄帶干擾的精準(zhǔn)頻率為

利用上述方法分別對窄帶干擾特征區(qū)域Ⅰ、Ⅱ和Ⅲ中窄帶干擾進行頻率估計，得到窄帶干擾頻率估計值如表3所示。從表3中可以看出，本文所提Candan算法可以精準(zhǔn)估計各窄帶干擾的頻率，說明本文方法可以有效解決傅里葉頻率估計算法中柵欄效應(yīng)的問題，同時該方法無需加窗處理和多次迭代計算，執(zhí)行效率較高。

表3 仿真PD信號的窄帶干擾頻率估計結(jié)果Table 3 Narrowband interference frequency estimation of simulated PD signal

1.4 快速獨立分量分析

染噪PD信號可以被視為PD信號和窄帶干擾的線性混合，進一步可以被視為PD信號、正弦信號和余弦信號的線性混合，因此本文采用快速獨立分量分析算法對染噪PD信號進行分離，該算法具體原理可參考文獻[20]。

由于快速獨立分量分析算法在開展混合信號的分離時，必須要求混合信號的數(shù)目大于源信號的數(shù)目，因此為了保證快速獨立分量分析算法可以有效分離出原始PD信號，需要利用1.3節(jié)頻率估計結(jié)果對混合信號矩陣進行擴展。對于每個窄帶干擾，構(gòu)建2個混合信號分別為sin(2π)和cos(2π)，和染噪PD信號共同構(gòu)成混合信號矩陣。

由于快速獨立分量分析算法自身存在的缺陷，所以分離的信號存在順序、幅值和極性的隨機性，即在分離的信號中難以確定哪組信號為PD信號，同時分離的PD信號和原始PD信號的幅值和極性可能存在差異。為了解決上述問題，本文在染噪PD信號前端添加長度為300的校正信號z(i)，該校正信號為

式中max(*)是求數(shù)組的最大值。

圖3 混合信號Fig. 3 Mixed signal

對于分離信號中PD信號的判斷問題，本文引入絕對相關(guān)數(shù)作為評判標(biāo)準(zhǔn)，定義各組分離信號前端長度為300的信號為r(i)，得到r(i)和z(i)的絕對相關(guān)數(shù)ρ為

取絕對相關(guān)數(shù)最大的一組分離信號作為分離后的PD信號(i)，然后利用式（11）對信號(i)進行修正，使得窄帶干擾抑制后的PD信號q(i)更接近于原始PD信號。

2 PD信號的窄帶干擾抑制

2.1 干擾抑制流程

本文所提PD信號中窄帶干擾抑制方法的具體步驟為：

1）對染噪PD信號進行廣義S變換，利用PD信號和窄帶干擾的時頻分布差異，在廣義S變換中劃分出窄帶干擾特征區(qū)域。

2）在各窄帶干擾的特征區(qū)域中，利用Candan算法對窄帶干擾頻率進行精準(zhǔn)估計。

3）利用估計得到的窄帶干擾頻率，構(gòu)建混合信號sin(2π)和cos(2π)。

4）在染噪PD信號的前端添加校正信號，并在混合信號sin(2π)和cos(2π)的前端添加零值信號，并將上述信號共同構(gòu)成混合信號矩陣。

5）在利用快速獨立分量分析算法對混合信號矩陣進行分析后，根據(jù)絕對相關(guān)數(shù)確定哪組信號為PD信號，同時利用校正信號修正分離的PD信號，得到窄帶干擾抑制結(jié)果。

2.2 干擾抑制效果

利用本文所提PD信號中窄帶干擾抑制方法對圖1（c）中染噪PD信號進行去噪，通過圖2的時頻分布圖可以確定窄帶干擾的數(shù)目為3，對應(yīng)的窄帶干擾特征矩陣也為3個，以此得到窄帶干擾抑制效果如圖4（a）所示。

為了說明本文方法的優(yōu)越性，使用奇異值分解方法[21]和FFT濾波方法[10]對圖1（c）中染噪PD信號進行窄帶干擾抑制，得到干擾抑制結(jié)果分別如圖4（b）和圖4（c）所示。從圖4中可以看出，本文所提方法可以有效抑制窄帶干擾，并且窄帶干擾抑制后的PD信號波形和原始PD信號波形的幅值和極性保持一致，說明本文方法可以有效恢復(fù)PD信號的波形特征，利于PD信號的后續(xù)分析。而傳統(tǒng)的奇異值分解方法和FFT濾波方法的窄帶干擾抑制結(jié)果中存在明顯的殘余噪聲，窄帶干擾抑制效果較差。

為了量化說明圖4中各窄帶干擾抑制方法的效果，本文引入3組參數(shù)對去噪效果進行評價，分別為信噪比系數(shù)SNR、波形相似系數(shù)NCC和均方誤差系數(shù)MSE[13]：

圖4 仿真PD信號的窄帶干擾抑制結(jié)果Fig. 4 Suppression results of narrowband interference for PD signal

計算得到圖4中各窄帶干擾抑制結(jié)果的評價參數(shù)值如表4所示。對比表4中各方法的評價參數(shù)可以看出，在3種方法中，本文所提方法能取得最高的SNR和NCC，并且能取得最低的MSE，說明本文所提方法可以很好恢復(fù)得到PD信號的波形，波形的畸變失真較少，利于后續(xù)的PD信號分析。

表4 窄帶干擾的抑制效果Table 4 Suppression effect of narrowband interference

需要說明的是，本文方法的計算復(fù)雜度較高，但是窄帶干擾的抑制效果更好，因此本文方法更加適合PD信號的離線分析。

3 實測信號的窄帶干擾抑制

為了分析本文所提方法的實際運用效果，在實驗室中制作缺陷電纜樣本，并開展工頻局部放電測試。電纜樣本選用10 kV交聯(lián)聚乙烯電纜，設(shè)置的缺陷為終端頭刀痕缺陷，具體的試驗接線圖如圖5所示。圖5中示波器的采樣率設(shè)置為100 MHz，采樣點數(shù)為1300。由于實驗室中電磁干擾較小，因此向采集得到的PD信號中疊加頻率分別為1.34、3.38和10.43 MHz，幅值分別為3、2和3 mV，相位分別為π/2、π/6和π/3 rad的窄帶干擾信號，得到實測的染噪PD信號如圖6所示。從圖6中可以看出，在窄帶干擾的影響下，PD信號已經(jīng)無法被識別提取，所以窄帶干擾的抑制是必要的。

圖5 局放測試接線示意圖Fig. 5 Wiring diagram of PD measurement

利用廣義S變換分析圖6中實測的染噪PD信號，得到染噪PD信號的時頻分析圖譜如圖7所示，從圖7中可以看出，窄帶干擾的數(shù)目為3，同時劃定各窄帶干擾對應(yīng)的時頻區(qū)域，得到窄帶干擾特征區(qū)域Ⅰ、Ⅱ和Ⅲ。在窄帶干擾特征區(qū)域內(nèi)，對窄帶干擾的頻率進行估計，得到實測染噪PD信號中各窄帶干擾的頻率如表5所示。

從表5中可以看出，本文所提方法可以精準(zhǔn)估計實測染噪PD信號的窄帶干擾頻率，利于后續(xù)窄帶干擾的抑制，繼續(xù)使用本文方法對窄帶干擾進行抑制，得到窄帶干擾抑制結(jié)果如圖8（a）所示。同時為了對比傳統(tǒng)的方法，進一步采用奇異值分解法和FFT濾波法對圖6中實測染噪PD信號進行窄帶干擾抑制，分別得到干擾抑制結(jié)果如圖8（b）和圖8（c）所示。對比圖8中各方法的窄帶干擾抑制結(jié)果可以看出：奇異值分解方法和FFT濾波方法的窄帶干擾抑制結(jié)果中存在明顯的殘余噪聲（圖8中圈內(nèi)標(biāo)記所示），窄帶干擾抑制效果有限；本文所提方法的窄帶干擾抑制效果明顯優(yōu)于傳統(tǒng)方法，能很好地抑制窄帶干擾。

表5 實測染噪PD信號的窄帶干擾頻率估計結(jié)果Table 5 Narrowband interference frequency estimation of measured PD signal with noise

圖6 實測的染噪PD信號Fig. 6 Measured PD signal with noise

為了進一步說明本文方法的窄帶干擾抑制效果，利用廣義S變換方法對本文方法的抑制結(jié)果進行時頻分析，得到對應(yīng)的時頻分布圖如圖9所示。對比圖7和圖9中可以看出，本文方法有效地抑制了窄帶干擾，對應(yīng)的處理結(jié)果中僅剩余白噪聲和PD信號，并且PD信號的時頻特征已經(jīng)清晰可見，可以進行相應(yīng)的提取分析。

圖7 實測染噪PD信號的時頻分布圖Fig. 7 Time-frequency distribution picture of measured PD signal with noise

圖9 窄帶干擾抑制結(jié)果的時頻分布Fig. 9 Time frequency distribution picture of narrowband interference suppression result

4 結(jié)論

1）廣義S變換可以對染噪PD信號進行時頻分析，根據(jù)PD信號和窄帶干擾的時頻特征，可以確定窄帶干擾的特征區(qū)域，從而對窄帶干擾特征量進行分析。

2）在各窄帶干擾的特征區(qū)域中，利用Candan算法可以準(zhǔn)確估計各窄帶干擾的頻率值，并且該方法無需多次迭代和加窗賦權(quán)，計算簡單。

3）利用帶校正信號的快速獨立分量分析方法對染噪PD信號中窄帶干擾進行抑制，可以解決分離信號的順序、幅值和極性的隨機性問題，能有效恢復(fù)PD信號的波形特征。

4）仿真測試結(jié)果中，本文方法的SNR=23.4377，NCC=0.9977，MSE=0.0031，說明本文所提方法可以很好恢復(fù)得到PD信號的波形，波形的畸變失真較少。實際測試的分析結(jié)果表明，相比于傳統(tǒng)的奇異值分解方法和FFT濾波法，本文方法的窄帶干擾抑制結(jié)果中殘余干擾能量較小，干擾抑制效果更明顯。