基于MATLAB軟件平臺的局放信號波形分析

蔣池劍 余祖良 鄭智勇

摘 要：局部放電是引起電力設備絕緣缺陷的主要因素，準確檢測局部放電信號對于提高電氣設備的可靠性至關重要。然而，局部放電信號中往往夾雜著大量噪聲，嚴重影響判斷。鑒于此，采用MATLAB軟件編程，利用希爾伯特-黃變換算法，通過主體函數EMD（y），實現對局部放電真實信號的提取。實例檢測結果表明，該方法具有一定的有效性和可靠性。

關鍵詞：希爾伯特-黃變換;MATLAB軟件;局部放電;信號提取

中圖分類號：TM93? 文獻標志碼：A? 文章編號：1671-0797（2022）07-0042-04

DOI：10.19514/j.cnki.cn32-1628/tm.2022.07.011

0? ? 引言

局部放電現象多出現于高壓電氣設備中，這種放電會造成導體間的絕緣局部短接，加速電氣設備絕緣介質的老化和損壞，從而縮短設備的使用壽命，嚴重時還會引發絕緣擊穿，損壞設備，造成人員生命、財產損失。然而，局放放電量小，發熱不明顯，很難對其識別和判斷。當局放現象明顯時，局放的發熱量往往已經造成一定的絕緣材料損失。

局部放電是電力變壓器中絕緣缺陷的主要特征之一，且由于變壓器在電力系統中的重要地位，準確檢測局部放電成為延長變壓器使用壽命、提高電力系統穩定性的重要措施。目前，國內外采用的局部放電識別方法主要有超聲波法、電氣法和超高頻法等。這些方法各有千秋，但正確應用的前提都是提取出局部放電的信號。然而，局部放電信號中往往夾雜著大量無關噪聲，嚴重影響判斷[1]。對此，1998年，中國臺灣海洋學家黃鍔在希爾伯特變換的基礎上提出了希爾伯特-黃變換[2]，依靠這種變換方法可以較為方便準確地提取和分離故障信號、干擾信號。黃鍔教授在完成對希爾伯特-黃變換的研究后，上傳了其變換的程序包，但程序包難以解碼，內部程序不可見，對于其他學者而言難以進行進一步的優化，故本文嘗試獨立于該程序包，利用MATLAB撰寫希爾伯特-黃內核變換代碼，并利用其對局放波形進行分析和研究。

1? ? 希爾伯特-黃變換算法

希爾伯特-黃變換（HHT）算法是由中國臺灣海洋學家黃鍔于1998年提出的，其基于希爾伯特變換并有所改進。該方法分為經驗模態分解（EMD）和希爾伯特譜分析（HSA）[3-4]，其處理信號的方式為：先對該信號進行經驗模態分解，分解出符合一定條件的模態因子（IMF），然后對所有模態因子進行希爾伯特譜分析，得到相應的希爾伯特譜，最終將所有的希爾伯特譜匯總，便能得到原始信號的希爾伯特譜。

HHT會對原始信號y（t）進行檢測，若其不單調，便對其進行經驗模態分解。其主要步驟為，首先繪出該信號的上下包絡線bs和bx，然后計算包絡線的均值j1：

j1=（1）

隨后用原信號減去均值，得到一個新的信號y11（t）。

y11（t）=y（t）-j1（2）

對得到的y11（t）進行檢測，檢測條件如下：（1）信號的極值點數和信號的零點數相差不超過1;（2）信號的上下包絡線的局部均值為0。若其不滿足以上條件，重復m次上述步驟，直到其滿足以上條件。得到滿足條件的信號y1m（t），其被稱為模態因子IMF。得到一個IMF后，令y1m（t）=r1（t）并記錄，然后用原信號減去該IMF，所得的殘余信號稱為c1（t）。若c1（t）仍不單調，則再次對其進行EMD分解，直到滿足IMF條件并得到新的ri（t）和殘余信號ci（t）。重復該步驟n次，直至分解出的cn（t）單調，并得到n個ri（t）和一個單調信號cn（t），該單調信號被稱為殘余信號。至此，原信號被分解為n個IMF和一個殘余信號，即：

y（t）=ri（t）+cn（t）（3）

再將這n+1個信號進行希爾伯特變化，得到瞬時頻率、瞬時能量、信號譜等，將所得的結果輸出，即完成信號的分解。

其核心為EMD處理。從頻域上來看，是高頻逐漸到低頻;從時域上來看，其特征時間尺度逐漸擴大。同時，在某些情況下，為了防止過度分解，往往需對其設置一個尺度，該尺度稱為SD，如式（4）所示：

SD=（4）

一旦計算出的SD小于某個閾值，則認為該次分解已經足夠小，并立即結束該次EMD計算。一般將閾值設置為0.2～0.3。

2? ? 希爾伯特-黃變換的代碼實現

了解了希爾伯特-黃變換的原理后，采用MATLAB編程實現這一變換。其主體函數如下EMD（y）所示，y是輸入信號。實現希爾伯特-黃變換的重點和難點在于如何實現EMD變換。

function a = EMD（y）

while ～dangdiao（y）

y1=y

sd=Inf

while （sd > 0.1） | ～imf（y1）

s1 = getspline（y1）;

s2 = -getspline（-y1）

h = y1-（s1+s2）/2

sd = sum（（y1-h）.^2）/sum（y1.^2）

y1=h

end

a{end+1} = y1

y=y-y1

end

end

在EMD（y）函數中包括3個副函數，函數名分別為“dangdiao” “imf” “getspline”。其中函數“dangdiao”用來判斷輸入信號是否為單調函數，若不為單調函數則進行循環;函數“imf”用來判斷信號是否滿足IMF條件，若不滿足則繼續循環;函數“getspline”為繪制上下包絡線函數。而sd作為限制尺度，用于防止過度分解。

函數“dangdiao”用以判斷信號在一定區間內是否為常數，其代碼如下：

function u=dangdiao（y）

u1 = jizhidian（y）;

if u1 > 0

u = 1;

else

u = 0;

end

end

函數“jizhidian”為極值檢測函數，能夠測量出原始信號的單調方式。若原始信號既有單調增區間，又有單調減區間，則返回值為0;除此之外，返回值為1。

function n=jizhidian（y）

k=find（diff（y）>0）

m=find（diff（y）<0）

k1=size（k）

m1=size（m）

if m1==0 | k1==0

n=1

else

n=0

end

end

函數“imf”函數的原理是測量原始信號的零點和極點的數量。若零點和極點的數量相差超過1，則不滿足IMF條件，返回值為1，否則返回值為0。其中，“peaks”函數作用為尋找信號全部的極值點，其將信號求二階導，若一階導值大于0，而二階導值小于0，則其必是極大值點。由此尋找出信號全部的極值點。

function u=imf（y）

N=length（y）

u1=sum（y（1：N-1）.*y（2：N） < 0）

u2=length（peaks（y））+length（peaks（-y））;

if abs（u1-u2）>1

u=0

else

u=1

end

end

函數“getspline”基于“spline”函數進行包絡線的繪制。“spline”是一個較為復雜的函數，其能從輸入數據中尋找內在聯系，并以此繪制出包絡線。

function s=getspline（y）

N1=length（y）;

p=peaks（y）;

s=spline（[0 p N1+1]，[0 y（p） 0]，1：N1）;

end

運行上述代碼，就可以實行信號分析，提取出原始信號的主體部分、次要信號以及噪聲信號。

3? ? 二維希爾伯特-黃變換的代碼實現

與一維信號EMD變換不同，BEMD變換主要面向二維圖片或波形，其處理流程與EMD變換類似，但增加了信號識別和信號特征提取，具體處理流程如圖1所示。

BEMD變換代碼實現函數主要由三部分組成：主體函數、包絡線繪制函數和單調判斷函數。其主體函數用于儲存分離出的IMF和余項;而包絡線繪制函數通過對函數進行描點處理，尋找其極值點，并以此繪制包絡線;同樣，單調判斷函數也通過對函數的描點處理，判斷其是否單調。其代碼如下：

clear all

img=imread（'wave.jpg'）

[r，c，d] = size（img）;

if d ～= 1

img = im2double（rgb2gray（img））;

else

img = im2double（img）;

end

imshow（img）

k=1;

m=0;

input_img=img;

A={}

while（～dangdiao（input_img））

[imf_de res_de]=baoluoxian（input_img）;

A{k}=imf_de;

input_img=res_de;

k=k+1

end

A{k}=res_de;

function? [imf_de res_de]=baoluoxian（input_img）

[width height]=size（input_img）;

x=1：width;

y=1：height;

input_img_temple=input_img;

SD=Inf

while（1）

[zmax imax zmin imin]=extrema2（input_img_temple）;

[xmax ymax]=ind2sub（size（input_img_temple），imax）;

[xmin ymin]=ind2sub（size（input_img_temple），imin）;

if SD<0.2

break

else

[zbs，～，～]=gridfit（ymax，xmax，zmax，y，x）;

[zbx，～，～]=gridfit（ymin，xmin，zmin，y，x）;

zbav=（zbs+zbx）./2

imf_de=input_img_temple-zbav;

SD=sum（sum（imf_de-input_img_temple）.^2）/sum（sum（imf_de）.^2）;

input_img_temple=imf_de;

end

res_de=input_img-imf_de;

end

end

function m=dangdiao（input_img）

[zmax imax zmin imin]=extrema2（input_img）;

[xmax ymax]=ind2sub（size（input_img），imax）;

[xmin ymin]=ind2sub（size（input_img），imin）;

xx=size（xmax，1）

xn=size（xmin，1）

u1=xx*xn

if u1>0

m=0

else

m=1

end

end

4? ? 實例測試

圖2為某變電站10 kV開關柜的局部放電信號波形圖，采用本文方法利用MATLAB對其進行信號提取和分解，得到的真實信號如圖3所示，次要信號如圖4所示，噪聲信號如圖5所示。本文方法能很好地將真實信號從原始信號中剝離出來。

5? ? 結語

本文研究了HHT的變換方式，利用MATLAB軟件代碼逐步實現了該變換，并取得了一定的成果。相比于簡單的希爾伯特變換，希爾伯特-黃變換增加了EMD處理，能夠更好地處理非線性復雜信號，且相較于常用的傅里葉變化和小波變換，其擁有較強的自適應性，泛用性更好。在出現奇異信號（不規則信號）時，主要采用積分分解法的傅里葉變化和小波變換會出現巨大的計算量，且容易產生虛假信號，嚴重影響判斷;而希爾伯特-黃變換能夠很好地處理奇異信號。

[參考文獻]

[1] 羅新，牛海清，宋廷漢，等.基于S變換和概率神經網絡的局部放電特征提取及放電識別方法[J].南方電網技術，2020，14（7）：17-23.

[2] 李光輝.基于希爾伯特黃變換及其改進方法的信號分析研究與應用[D].成都：成都理工大學，2012.

[3] 曾祥，周曉軍，楊辰龍，等.基于經驗模態分解和S變換的缺陷超聲回波檢測方法[J].農業機械學報，2016，47（11）：414-420.

[4] 鄔蒙蒙，周懷來，林萍，等.改進的完備經驗模態分解和廣義S變換相結合的地震信號衰減分析[J].地球物理學進展，2020，35（5）：2001-2008.

收稿日期：2022-01-11

作者簡介：蔣池劍（1978—），男，浙江人，高級工程師，研究方向：抽水蓄能電站生產管理。

通信作者：鄭智勇（1978—），男，湖南人，工程師，從事電氣設備管理工作。