李 晨，李 川，姜 飛，張長勝

(昆明理工大學(xué) 信息工程與自動化學(xué)院，云南 昆明 650500)

隨著電力系統(tǒng)的發(fā)展和廣泛應(yīng)用，電能對人們的生產(chǎn)和生活都有不可替代的作用。但是在電力系統(tǒng)中，由于大量電子設(shè)備和非線性負載的存在，使相關(guān)的電流和電壓波形產(chǎn)生較大程度的畸變，即諧波污染[1-3]。近年來，諧波對電網(wǎng)的影響也越來越受到電力部門的重視，國內(nèi)外的專家、學(xué)者均對電力諧波的影響[4-5]、計量[6-7]和抑制[8]進行廣泛的研究。1999年，張伏生等人[9]，通過FFT(Fast Fourier Transform，F(xiàn)FT)加窗插值方法提高諧波測量的精度，但都存在一定程度上的頻譜泄露和柵欄效應(yīng)。2012年，房國志等人[10-11]，提出基于FFT和小波包變換的綜合檢測方法，利用傅里葉變換識別諧波分量，再使用小波包對諧波進行定位和提取，有效的減少了信號提取的運算量，但對于高頻信號分辨率較低易造成信號頻帶混疊現(xiàn)象。

本文針對HHT在進行信號處理時所出現(xiàn)的端點效應(yīng)問題和模態(tài)混疊問題提出了改進方法，采用包絡(luò)極值延拓法來減少端點“飛冀現(xiàn)象”，通過對原始信號加入白噪聲輔助分析來消除模態(tài)混疊現(xiàn)象。最后，進行實驗仿真。結(jié)果表明，改進的希爾伯特-黃變換方法對于暫態(tài)擾動信號有更精確的分析和定位。

1 HHT變換原理

希爾伯特-黃變換(Hilbert-Huang Transform，HHT)是1998年Norden E.Huang等人[12]提出的一種經(jīng)驗數(shù)據(jù)分析方法，并引入Hilbert譜的概念和Hilbert譜的分析方法。HHT是分析非線性非平穩(wěn)信號的時頻分析方法，其中包括兩個部分。第一部分為經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)假設(shè)任何復(fù)雜信號都可以分解為有限數(shù)目且具有一定物理定義的固有模態(tài)函數(shù)(Intrinsic Mode Function，IMF)，并且能夠依據(jù)信號特點將信號從高到低分解成不同頻率的一系列IMF，在自適應(yīng)獲取信號基函數(shù)的同時也造成了模態(tài)混疊問題。第二部分為Hilbert譜分析(Hilbert Spectrum Analysis，HSA)，將用Hilbert求得的每一階IMF的瞬時頻率和瞬時幅值，既Hilbert譜。IMF需滿足兩個條件[13-14]：(1) 信號極值點的數(shù)量與零點數(shù)相等或相差是1；(2) 信號的由極大值定義的上包絡(luò)和由極小值定義的下包絡(luò)的局部均值為零，即IMF的上下包絡(luò)線關(guān)于時間軸對稱。

1.1 經(jīng)驗?zāi)B(tài)分解(EMD)

(1) 給定任意信號f(t)，得到信號的極值點(至少兩個)，用三次樣條插值法分別擬合信號f(t)的極小值和極大值得到上下兩條包絡(luò)線v(t)和u(t)，計算上下兩條包絡(luò)線每點對應(yīng)的均值mi(t)即mi(t)=[ui(t)+vi(t)]/2；

(2)在每點所對應(yīng)的fi(t)，用fi(t)減去每點所對應(yīng)的均值mi(t)得到hi(t)=fi(t)-mi(t)。

用得到的hi(t)檢驗是否滿足IMF的兩個條件，若滿足則hi(t)為f(t)的第一個IMF即IMF1；若不滿足則hi(t)作為原始數(shù)據(jù)fi(t)重復(fù)步驟(1)和(2)得到hij(t)，直到hij(t)滿足IMF的兩個條件。

(3)將信號f(t)用EMD分解得到的第一個IMF分量c1，r(t)作為新的原始數(shù)據(jù)得到r(t)=f(t)-c1。

重復(fù)步驟(1)、(2)、(3)得到信號f(t)的第2個IMF分量c2，循環(huán)n次得到信號f(t)的第n個IMF分量

(1)

式中rn代表函數(shù)的平均趨勢，是信號的殘余分量，其曲線為單調(diào)函數(shù)或定值。各個IMF分量cn表示了信號的特征尺度，體現(xiàn)非線性非平穩(wěn)信號的內(nèi)在模態(tài)特征。

1.2 Hilbert譜分析(HSA)

將式(1)中信號的IMF分量ci做希爾伯特譜變換

(2)

得到信號f(t)的解析信號Yi(t)=ci(t)+jHi(t)，由此得到瞬時頻率、瞬時相位和瞬時振幅為

(3)

(4)

(5)

所以f(t)的解析信號可表示為

Yi(t)=ci(t)+jHi(t)=ai(t)ejφi(t)=ai(t)ej?φidt

(6)

以上公式表示了信號f(t)頻率、相位和幅值之間的關(guān)系。信號f(t)的幅值可以用時間和頻率的函數(shù)H(ω,t)表示

(7)

對時間進行積分即可得到信號的邊際譜

(8)

其中，H(ω,t)是對整個信號的幅值在頻率上隨時間變化規(guī)律的描述；而h(ω)是信號在頻率上的總振幅(總能量)。

2 HHT在電力諧波中的應(yīng)用

基于實際電網(wǎng)中諧波源的不斷變化，造成諧波信號的次數(shù)、幅值也隨之變化。本文采用高次衰減電壓突變諧波信號模型f(t)[15-16]。

f(t)=

(9)

其中，設(shè)采樣頻率為3 200 Hz，采樣點數(shù)N=2 048。在采樣周期內(nèi)，3次諧波在0.06 s處開始衰減幅值由217 V降到87 V，在0.14 s的時候恢復(fù)到217 V；9次諧波在0.16 s處開始衰減幅值由93 V降到56 V，在0.24 s恢復(fù)到93 V。原始信號波形及其頻譜圖如圖1所示。

圖1 原始信號波形及頻譜圖

由圖1信號的頻譜圖可知，原始信號中分別包含頻率為50 Hz、150 Hz、和450 Hz的基波、3次諧波和9次諧波。

2.1 模態(tài)混疊

對信號進行EMD實際就是把信號的頻率從高到低的一個分解過程。在分解過程中由于基波幅值分量遠大于諧波幅值分量，可能會產(chǎn)生尺度交叉現(xiàn)象，其結(jié)果會導(dǎo)致模態(tài)混疊現(xiàn)象。如圖2所示。

圖2 原始信號的EMD分解

由圖2可以看出，對原始信號進行EMD分解，生成各個IMF分量，其中IMF1的波形在0.16 s時產(chǎn)生了壓降在0.24 s時恢復(fù)，IMF2波形在0.06 s時，9次諧波衰減至3次諧波，在0.16 s時恢復(fù)。IMF3和 IMF4中兩端波形都有變形，在下面IMF的幅頻曲線圖中更直觀的表現(xiàn)出。

2.2 端點效應(yīng)

對原始信號進行EMD分解后所得到的IMF分量進行Hilbert譜分析時，無法保證端點恰好為極值點。如果在進行Hilbert譜分析時把端點當(dāng)作極值點處理，則在兩端會造成信號的失真，出現(xiàn)明顯的“飛冀”現(xiàn)象?？稍贗MF的瞬時頻譜圖中清晰的看出，如圖3所示。

圖3 各個IMF的瞬時頻率

3 改進HHT的仿真分析

3.1 包絡(luò)極值延拓法

通過包絡(luò)極值延拓法來消除分解過程中出現(xiàn)的端點“飛冀”現(xiàn)象。在分解IMF分量時，信號的極值點分別構(gòu)成上包絡(luò)線和下包絡(luò)線。因此為了消除端點現(xiàn)象，依據(jù)極值點包絡(luò)線的特征分別向左右兩端進行延拓，同時考慮IMF所需要滿足的兩個條件，所以分別向左右兩端延拓兩個極值點。

設(shè)有一組信號數(shù)據(jù)X(T)，其中X(T)中有m個極大值點，有n個極小值點，采樣周期為Ts，采樣點數(shù)為N。

左延拓：設(shè)左端的第一個包絡(luò)波形特征值為l1

(10)

得向左延拓極值的位置和數(shù)值分別為

Tm0=Tm1-l1T;Xm(0)=Xm(1)

(11)

Tm-1=Tm1-2l1Ts;Xm(-1)=Xn(1)

(12)

Tn0=Tn1-l1Ts;Xn(0)=Xn(1)

(13)

Tm0=Tm1-2l1Ts;Xn(-1)=Xn(1)

(14)

右延拓：設(shè)右端的第一個包絡(luò)波形特征值為l2

(15)

得向右延拓極值的位置和數(shù)值分別為

Tm+1=Tm+l2Ts;Xm(m+1)=Xm(m)

(16)

Tm+2=Tm+2l2Ts;Xm(m+2)=Xm(m)

(17)

Tn+1=Tn+l2Ts;Xn(n+1)=Xn(n)

(18)

Tn+2=Tm+2l2Ts;Xn(n+2)=Xn(n)

(19)

當(dāng)左右兩端的函數(shù)值小于鄰近端點的極值時，為防止極值在包絡(luò)線以外，需對其進行特殊處理

Tm0=T1,Xm(0)=X(1);X(1)>Xm(1)

(20)

Tn0=T1,Xn(0)=X(1);X(1)

(21)

Tm+1=Tn,Xm(n+1)=X(n);X(n)>Xm(m)

(22)

Tn+1=Tn,Xn(n+1)=X(n);X(n)

(23)

由于分別向左右延拓了4個極值點，最后在進行EMD分解時可以左右兩邊各去掉幾個采樣點數(shù)以消除延拓對實際信號的影響。

3.2 平均經(jīng)驗?zāi)B(tài)法

分解時產(chǎn)生的模態(tài)混疊現(xiàn)象一方面是和EMD算法有關(guān)，另一方面也和原始信號的頻率特征有關(guān)。當(dāng)原始信號的時間尺度存在階躍性變化時，則分解的IMF分量包含不同時間尺度特征分量。針對這一現(xiàn)象，本文采用平均經(jīng)驗?zāi)B(tài)法來抑制模態(tài)混疊現(xiàn)象。在進行EMD分解時相當(dāng)于一個二階濾波器組，當(dāng)信號出現(xiàn)間斷點等情況時濾波器無法對信號準(zhǔn)確分解。因此采用EEMD方法對信號進行噪聲輔助分析，噪聲輔助信號方法就是在原始信號中加入白噪聲來平滑脈沖干擾，而加入的白噪聲對各個固有模態(tài)分量的分解并沒有影響。白噪聲的頻譜均勻散落在信號整個時頻空間上，由于EEMD能夠自適應(yīng)對信號進行分解，所以不同時間特征的信號會自適應(yīng)的分布到合適的時間尺度上。并且由于加入的白噪聲具有零均值噪聲的特性，所以各個固有模態(tài)分量在經(jīng)過多次平均后可消除白噪聲。

3.3 改進HHT的仿真

采用改進的HHT算法對式(9)信號模型仿真進行對比。圖4為原始信號中加入均值為0、方差為1時正態(tài)分布的隨機數(shù)。

圖4 原始信號加噪聲

對加入隨機白噪聲的原始信號進行EEMD分解得到的各個IMF分量，由圖5和圖6可以看出加入的隨機白噪聲被很好的分解出來，原始信號中的高次衰減電壓突變也得到了很好的處理，IMF2是加入的噪聲的余量，其中IMF8的波形接近正弦波形。

圖5 信號的IMF1-IMF4分量

圖6 信號的IMF5～IMF8分量

圖7 信號的IMF1～IMF4的瞬時頻率

圖8 信號IMF5～IMF8的瞬時頻率

由圖7和圖8可以觀察到，在曲線的兩端端點“飛冀”問題，可通過包絡(luò)極值延拓法得到良好的抑制。

4 結(jié)束語

本文針對希爾伯特-黃算法在電能諧波分析中出現(xiàn)的端點效應(yīng)和模態(tài)混疊問題，提出用包絡(luò)極值延拓法和平均經(jīng)驗?zāi)B(tài)法來改善。利用Matlab建立一個高次衰減電壓突變的典型諧波信號模型，對該模型進行仿真。通過改進前后的希爾伯特-黃算法，比較經(jīng)驗?zāi)B(tài)分解得到的各個IMF分量及其幅頻曲線，改進后的希爾伯特-黃算法明顯改善了端點效應(yīng)和模態(tài)混疊問題。