楊 彥 ,王 浩 ,曹德健,趙 力*
(1.鹽城紡織職業(yè)技術(shù)學(xué)院,江蘇 鹽城 650500;2.東南大學(xué)信息科學(xué)與工程學(xué)院,南京 210096)
由于電力系統(tǒng)的復(fù)雜性,在遠(yuǎn)程傳輸以及運(yùn)載中會產(chǎn)生的不平穩(wěn)的故障信號,這也就使得信號突變的檢測在電力系統(tǒng)中占有重要的地位。其常用的傳統(tǒng)數(shù)學(xué)工具為短時傅里葉變換,短時傅里葉變換將信號劃分為許多小的時間間隔,用傅里葉變化分析每個時間間隔,一邊確定該時間間隔的頻率。然而對于短時傅里葉變換,時間—頻率窗口的寬度對于所有頻率的譜具有不變特性,這一點(diǎn)不適合于非平穩(wěn)信號的高頻與低頻部分的特性分析。后來很多工程應(yīng)用者嘗試將小波變換理論引入電力系統(tǒng)故障分析中,因為小波變換能夠在時域和頻域同時得到比較高的分辨率,然而對信號進(jìn)行小波變換必須預(yù)先確定小波的基函數(shù),且信號分解的效果取決于小波基函數(shù)的選擇,并不能夠保證最優(yōu)的小波分解效果[1]。1998年,美國國家宇航局(NASA)的Norden和E.Huang 等人在對瞬時頻率的概念進(jìn)行深入研究的基礎(chǔ)上,創(chuàng)立了希爾伯特-黃變換,該方法被認(rèn)為是近年來對已傅里葉變換為基礎(chǔ)的線性和穩(wěn)態(tài)譜分析的一個重大突破[1]。它是一種基于經(jīng)驗?zāi)J椒纸釫MD(Empirical Mode Decomposition)的時頻分析處理方法。它將非平穩(wěn)信號按不同尺度的波動或趨勢逐級分解成若干個本征模式分量IMF(Intrinsic Mode Function),然后對IMF 進(jìn)行Hilbert 變換,得到信號的時頻譜圖具有良好的時頻分辨率,從而具有良好的局部視頻特性表現(xiàn)能力,因此能夠更加準(zhǔn)確地反映原信號的特性,可以用來檢測突變信號。該方法吸取了小波變換的多分辨率的優(yōu)勢,同時克服了小波變換中選取小波基的困難。能夠更加準(zhǔn)確地檢測出突變信號,從而正確定位出信號的故障點(diǎn)。
電力系統(tǒng)的不斷發(fā)展,使得輸送電壓和容量的逐步提高,有輸電線故障造成的損失也越來越大。電力系統(tǒng)主要有短路故障、斷路故障、復(fù)雜故障。用Simulink 中的Power system 仿真工具包來進(jìn)行故障信號仿真,通過仿真可以得到如圖1 幾種故障的A相電壓信號波形圖。
圖1 幾種故障的電壓波形圖
經(jīng)驗?zāi)J椒纸馐?998年Huang 提出來的一種全新的信號處理方法[2-3]。它通過對信號的篩選將信號分解成不同頻率的IMF,IMF 具有以下特點(diǎn):第1 IMF 的極值(包含極大值和極小值)數(shù)與過零點(diǎn)的數(shù)目相等或最多相差一個;第2 在任意時刻,IMF 的上包絡(luò)線與下包絡(luò)線的均值必須是零。具體的分解過程如下[4]:
第1 步:首先找到待分析信號x(t)所有的極大值點(diǎn)并將其使用3 次樣條函數(shù)插值(在下面做了描述)成為原數(shù)據(jù)序列的上包絡(luò)線;找到x(t)所有的極小值點(diǎn)將其使用3 次樣條函數(shù)插值成為原數(shù)據(jù)序列的下包絡(luò)線,然后求出上包絡(luò)線和下包絡(luò)線的平均包絡(luò)線m1(t);
第2 步:將原數(shù)據(jù)序列減去平均包絡(luò)后即可得到一個去掉低頻的新數(shù)據(jù)序列h1(t)=x(t)-m1(t);然后判斷h1(t)是否滿足IMF 的條件,若不滿足,將h1(t)看作新的x(t),重復(fù)上述處理過程,直到h1(t)滿足IMF 條件時,記c1(t)=h1(t),視為IMF1;
第3 步:將r(t)=x(t)-h(huán)1(t)看作新的x(t),重復(fù)以上(1)和(2)的步驟,即可以依次得到IMF2,IMF3,…,直到cn(t)或r(t)滿足給定的終止條件時篩選結(jié)束。
按照上面的分解步驟,原始數(shù)據(jù)序列x(t)可以由這些IMF 分量及一個均值或者趨勢項r(t)來表示:
其中每一個IMF 分量都代表一個特征尺度的數(shù)據(jù)序列,所以可以說EMD 過程實際上是將原數(shù)據(jù)序列分解成為不同特征波動的序列。
其中Re 表示取實部,ai(t)表示每一個本征模函數(shù)經(jīng)希爾伯特變換后的瞬時振幅,fi(t)表示每一個本征模函數(shù)經(jīng)希爾伯特變換后的瞬時頻率。若以時間、瞬時頻率為自變量,幅度以等高線的形式表示,則可以將三者之間的關(guān)系表示為一個三維時頻譜,又稱為Hilbert 幅度譜,簡稱為Hilbert 譜,記為H(f,t),則有
下面進(jìn)一步來分析Hilbert 譜,定義邊緣譜為h(f),有
其中T 表示信號序列的時間長度。根據(jù)邊緣譜的表達(dá)式,可以看出,邊緣譜是對三維時頻譜在時間上進(jìn)行積分,便形成了只有頻率和幅值的二維的關(guān)系圖,邊緣譜表征了整組數(shù)據(jù)在每個頻率點(diǎn)的積累的幅度值分布,而傅里葉頻譜在每一點(diǎn)的幅度值表征了信號里有一個含有此頻率的三角函數(shù)的組分,傅里葉頻譜和邊緣譜是不同意義的兩種頻譜與幅值關(guān)系表達(dá)式。
從上述理論,我們可以看出,EMD 算法是具有自適應(yīng)性的,主要表現(xiàn)在:生成的基函數(shù)具有自適應(yīng)性,IMF 的頻率分辨率也是具有自適應(yīng)性的,且同時具有濾波特性[5]。所以可以把EMD 分解看做是一組具有自適應(yīng)特性的帶通濾波器,其截至頻率和帶寬隨被分解信號的不同而不同[6-7]。同小波變換相比,EMD 是有自適應(yīng)的,而小波分析則沒有,因為小波的分解尺度被選定,分解得到的時域波形就是固定頻率段的,且頻率段只與分析頻率有關(guān),與信號本身無關(guān)[8]。這就是EMD 相比于小波變換的優(yōu)點(diǎn)。
EMD 分解的一個重要步驟是構(gòu)造信號的上下包絡(luò)線,從而得到信號的瞬時平均包絡(luò)。Huang 等提出的以信號的極大值點(diǎn)來擬合上包絡(luò)線,極小值點(diǎn)來擬合下包絡(luò)線,其中擬合算法采用的是3 次樣條插值方法。然而,我們知道,在信號的開始端和結(jié)束端往往不是極大值或者極小值點(diǎn),因此3 次樣條曲線容易在數(shù)據(jù)的兩個端點(diǎn)出現(xiàn)發(fā)散現(xiàn)象,且EMD的過程也使得這種發(fā)散現(xiàn)象不只局限于端點(diǎn)處的小范圍,而是伴隨著每次的篩選過程逐漸往內(nèi)傳播,會引起最終結(jié)果較大的失真,這種現(xiàn)象稱為“端點(diǎn)效應(yīng)”。此外,我們知道,EMD 的過程就是對信號進(jìn)行篩選的過程,篩選對信號有兩個結(jié)果:第1,駝峰減少,使得極大值為正數(shù)或零,極小值為負(fù)數(shù)或零;第2,光滑了非對稱幅度。然而過度的篩選會使得信號在幅度上失去物理意義,且當(dāng)前的電力設(shè)備很多都不能滿足這種過度的篩選迭代,在這里將這種現(xiàn)象稱為過度篩選現(xiàn)象。
在電力系統(tǒng)的復(fù)雜環(huán)境下,信號波形可能含有多個諧波和噪聲,端點(diǎn)效應(yīng)影響程度可能會加大,所以需要采取一定的方法來改善這種端點(diǎn)效應(yīng)造成的誤差。目前抑制這種效應(yīng)主要采用兩種方法:第1,尋找更加精確的求解上下包絡(luò)的方法來替代已有的3 次樣條插值方法,比如有約束3 次樣條插值法,分段冪函數(shù)法,高次樣條插值法;第2,在信號的兩端選取合適的極大值和極小值點(diǎn),使上下包絡(luò)線能夠完整的包絡(luò)整個信號。在這里,我們采用第2種方法來減少端點(diǎn)效應(yīng)所帶來的誤差,由于電力信號在一般情況下非故障點(diǎn)表現(xiàn)為一系列正弦信號之和,故電力信號在非故障處可以看作周期信號,本文將端點(diǎn)處的值作為極大值或極小值進(jìn)行處理,當(dāng)端點(diǎn)值為正數(shù)時,將其作為極大值,對其取反作為極小值;當(dāng)端點(diǎn)值為負(fù)數(shù)時,將其作為極小值,對其取反作為極大值。采用這種方法信號的包絡(luò)線得到了很大的改善,能夠減少端點(diǎn)效應(yīng)所帶來的誤差,能夠滿足大多數(shù)電力設(shè)備的要求。
對于過度篩選現(xiàn)象,本文降低了篩選的要求,即降低了IMF 的要求。我們知道,當(dāng)一個信號完全滿足IMF 的要求時,其上下包絡(luò)線是關(guān)于時間軸完全對稱的,包絡(luò)線均值是處處為零的。一般情況,即使在誤差范圍內(nèi),也需要進(jìn)行很多次的篩選,在這里,Huang等人使用esd標(biāo)準(zhǔn)來降低篩選條件。esd的表達(dá)式為
其中hi-1(t)和hi(t)分別代表經(jīng)過第i-1 次和第i次篩選后的信號。一般情況下,當(dāng)esd下降到0.2或者0.3 時,結(jié)束篩選過程。然而在esd中分母h2i-1(t)可能會出現(xiàn)零的情況,一般情況下這種為零現(xiàn)象是不可避免的,因此,esd的計算可能是不準(zhǔn)確的。
從上面所述,當(dāng)信號滿足IMF 的要求時,其上下包絡(luò)線是關(guān)于時間軸完全對稱的,所以在篩選的過程中,信號的上下包絡(luò)線是關(guān)于時間軸趨向?qū)ΨQ的,那么我們可以根據(jù)包絡(luò)線的對稱程度來決定是否結(jié)束篩選。在這里,我們設(shè)第i 次的上下包絡(luò)線數(shù)據(jù)序列分別為mup(n)和mdown(n),為了表示兩個數(shù)據(jù)變量的近似程度,我們采用了相關(guān)系數(shù)標(biāo)準(zhǔn)。根據(jù)相關(guān)系數(shù)的定義,設(shè)ri為第i 次時的上下包絡(luò)線相關(guān)系數(shù),表達(dá)式為
在一般情況下,我們將相關(guān)系數(shù)ri設(shè)置為大于等于0.985 時結(jié)束篩選。這樣可以有效降低EMD算法的循環(huán)次數(shù)。
圖2為三相短路時的電壓信號波形,我們在0.03 s 時設(shè)置故障,在0.065 s 時切除。
圖2 三相短路時的電壓信號波形
應(yīng)用前兩節(jié)方法對該信號進(jìn)行EMD 分解,使用Matlab 仿真,結(jié)果為圖3。
圖3 信號的EMD 分解
對IMF 做希爾伯特變換,可得到圖4。
圖4 希爾伯特變換后的時頻關(guān)系圖
通過在Matlab 中的放大觀察,我們看到信號在0.030161和0.065125 出的頻率很高,也就是故障發(fā)生和切除的時間,其誤差0.125 ms~0.161 ms,故障時間定為準(zhǔn)確,且運(yùn)行時間在許可的時間內(nèi)。
電力系統(tǒng)對于繼電保護(hù)裝置的要求是快速,準(zhǔn)確。本文采用了EMD 方法將信號分解,根據(jù)其IMF分量的希爾伯特變換得到信號的時頻特性,從而能夠準(zhǔn)確定位故障發(fā)生的時間,在下一步的工作中,我們應(yīng)該尋求對故障信號進(jìn)行處理分類,使得繼電保護(hù)器能夠迅速的分辨出故障的類型,以便做出相應(yīng)的動作。
[1]李天云,趙妍,李楠.基于EMD 的Hilbert 變換應(yīng)用于暫態(tài)信號分析[J].電力系統(tǒng)自動化,2005,29(4):49-52.
[2]Huang N E,Shen Z,Long S R.The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Nonstationary Time Series Analysis[C]//Proc.R.Soc.London A,1998,4:903-995.
[3]Huang N E,Stever R L.A New View of Nonlinear Water Waves:the Hilbert Spectrum[C]//Annual Review of Fluid Mechanics.1999,1:417-457.
[4]Balocchi R,Menicucci D,Varanini M.Empirical Mode Decomposition to a Approach the Problem of Detecting Sources from a Reduced Number of Mixtures[C]//Cancun,Mexico:Proceedings of the 25th Annual International Conference of the IEEE EMBS,2003.
[5]Raia V K,Mohanty A R.Bearing Fault Diagnosis Using FFT of Intrinsic Mode Functions in Hilbert-Huang Transform[J].Mechanical Systems and Signal Processing,2007,21(6):2607-2615.
[6]Wu Z H,Huang N E.Ensemble Empirical Mode Decomposition:A Nosie Assisted Data Analysis Method[R].Calverton Center for Ocean-Land-Atmosphere Studies,2005:855-895.
[7]Gabriel R,Patrick F.One or Two Frequencies? The Empirical Mode Decomposition Answers[J].IEEE Transactions on Signal Processing,2008,56(1):85-96.
[8]Zheng G T,McFadden P D.A Time-Frequency Distribution for Analysis of Signals with Transient Components and Its Application to Vibration Analysis Transaction of the ASME[C]//Proc.R.Soc.London A,July 1999,2:328-333.