汪頌軍，劉滌塵，廖清芬，周雨田，王亞俊，王乙斐，趙一婕

（武漢大學(xué) 電氣工程學(xué)院，湖北 武漢 430072）

0 引言

隨著全國(guó)聯(lián)網(wǎng)程度的不斷加深，低頻振蕩日益成為危及電網(wǎng)安全穩(wěn)定運(yùn)行的突出問(wèn)題之一?；谙到y(tǒng)日常運(yùn)行工況下的實(shí)測(cè)軌跡進(jìn)行的低頻振蕩模式在線辨識(shí)，對(duì)于低頻振蕩的在線監(jiān)測(cè)、預(yù)警、控制器設(shè)計(jì)、系統(tǒng)阻尼優(yōu)化等具有重要意義[1-3]。

目前對(duì)電力系統(tǒng)低頻振蕩進(jìn)行分析一般基于系統(tǒng)擾動(dòng)后的響應(yīng)軌跡，常用的方法有傅里葉變換、ARMA算法、Prony分析、小波算法、希爾伯特-黃變換（HHT）、旋轉(zhuǎn)不變技術(shù)（ESPRIT）等[4-6]，已有研究證明這些方法均能較準(zhǔn)確地辨識(shí)出系統(tǒng)的模式信息。但是這類(lèi)辨識(shí)方法均需對(duì)系統(tǒng)施加激勵(lì)，只能在系統(tǒng)發(fā)生振蕩后辨識(shí)出當(dāng)時(shí)的系統(tǒng)特性，作出告警，并不能給出系統(tǒng)正常運(yùn)行狀態(tài)下的模式信息，也無(wú)法實(shí)現(xiàn)“全天候”監(jiān)控和真正意義上的預(yù)警[7-8]。近年來(lái)基于環(huán)境激勵(lì)的模式辨識(shí)研究已在橋梁、建筑、汽輪機(jī)、飛機(jī)等[9]領(lǐng)域展開(kāi)，電力領(lǐng)域也逐步開(kāi)展了一些研究：文獻(xiàn)[10]探討了工況模式分析在低頻振蕩辨識(shí)中的應(yīng)用，論述了電力系統(tǒng)低頻振蕩和一般振動(dòng)力學(xué)數(shù)學(xué)模型的相似性；文獻(xiàn)[7]通過(guò)系統(tǒng)建模證明了負(fù)荷波動(dòng)作為激勵(lì)信號(hào)的可行性，以上研究為工況模式分析在電力系統(tǒng)的應(yīng)用做好了鋪墊。自然激勵(lì)技術(shù) NExT（Natural Excitation Technique）[9，11-12]是工況模式分析的有效方法，它可以避免傳統(tǒng)模式辨識(shí)人工激勵(lì)給電力系統(tǒng)帶來(lái)的傷害，實(shí)現(xiàn)系統(tǒng)在工況下的在線模式識(shí)別，其在電力領(lǐng)域僅在文獻(xiàn)[13]中有應(yīng)用，值得更進(jìn)一步的研究。

集合經(jīng)驗(yàn)?zāi)Ｊ椒纸釫EMD（Ensemble Empirical Mode Decomposition）[14-15]是經(jīng)驗(yàn)?zāi)Ｊ椒纸?EMD（Empirical Mode Decomposition）的發(fā)展，既繼承了EMD的優(yōu)點(diǎn)，同時(shí)通過(guò)向原始信號(hào)添加白噪聲改善了模式混疊問(wèn)題。它是一種自適應(yīng)的分析方法，適合處理工程中的非線性、非平穩(wěn)信號(hào)。

基于此，本文將NExT與EEMD結(jié)合，提出一種適合于工況下的低頻振蕩主導(dǎo)模式辨識(shí)方法。首先通過(guò)EEMD將工況下的實(shí)測(cè)時(shí)變信號(hào)分解為若干個(gè)單自由度的本征模式，憑借EEMD時(shí)空濾波器、互相關(guān)系數(shù)和能量權(quán)重比篩選出真實(shí)主導(dǎo)模式分量，然后通過(guò)NExT求互相關(guān)函數(shù)以代替?zhèn)鹘y(tǒng)的脈沖響應(yīng)，至此就實(shí)現(xiàn)了工況下把非平穩(wěn)的多自由度系統(tǒng)的模式辨識(shí)問(wèn)題轉(zhuǎn)化為單自由度系統(tǒng)的脈沖響應(yīng)辨識(shí)問(wèn)題，最后對(duì)互相關(guān)函數(shù)分別通過(guò)Teager能量算子TEO（Teager Energy Operator）求出模式幅值和頻率、能量分析求阻尼比、峰值法求相位。

1 EEMD

1.1 EEMD原理

為了改善EMD出現(xiàn)的模式混疊現(xiàn)象，Wu Zhaohua等提出了EEMD方法，其本質(zhì)是一種疊加高斯白噪聲的多次EMD。下面介紹EEMD的步驟。

a.產(chǎn)生N條隨機(jī)正態(tài)分布的白噪聲ni（t）分別加入到原始信號(hào) x（t）中，得每次加噪后的信號(hào) xi（t）為：

b.對(duì)每條加噪后的信號(hào) xi（t）進(jìn)行 EMD，得：

其中，cij（t）為第i次加入高斯白噪聲后分解得到的第 j個(gè)本征模式函數(shù)（IMF）分量；ri（t）為對(duì) xi（t）進(jìn)行EMD后的余項(xiàng)。

c.將N組分解結(jié)果中對(duì)應(yīng)的IMF分量求均值，則：

其中，cj（t）為對(duì)原始信號(hào)進(jìn)行EEMD得到的第 j個(gè)IMF 分量；r（t）為余項(xiàng)。

EEMD算法流程如圖1所示。

圖1 EEMD算法流程圖Fig.1 Flowchart of EEMD algorithm

1.2 EEMD時(shí)空濾波器

EEMD得到的IMF分量頻率由大到小，利用這點(diǎn)本文基于EEMD構(gòu)造時(shí)空濾波器組（TFB），其濾波信號(hào)為：

其中，x?（t）為濾波器輸出信號(hào)；l，h?[1，M]。

當(dāng)l=1且h＜M時(shí)，TFB為高通濾波器；當(dāng) l＞1且h=M 時(shí)，TFB為低通濾波器；當(dāng) 1＜l≤h＜M時(shí)，TFB為帶通濾波器。

1.3 互相關(guān)系數(shù)

原始信號(hào)經(jīng)EEMD所得的IMF分量中常包含一些虛假分量，而這些虛假分量與原始信號(hào)無(wú)關(guān)且不能反映原始信號(hào)的特征，因此有必要將其辨別出來(lái)并予以剔除?；诨ハ嚓P(guān)系數(shù)來(lái)判斷IMF分量的真?zhèn)问且环N行之有效的方法，通過(guò)求取分解后的各IMF分量與原始信號(hào)的互相關(guān)系數(shù)，假定互相關(guān)系數(shù)很小的對(duì)應(yīng)IMF分量為偽分量。設(shè)2個(gè)時(shí)間序列 x（n）、y（n），它們的互相關(guān)系數(shù) ρxy的表達(dá)式如式（6）所示。

互相關(guān)系數(shù)只是一個(gè)比率，不是等單位量度，互相關(guān)系數(shù)的正負(fù)號(hào)只表示相關(guān)的方向，絕對(duì)值表示相關(guān)的程度。

1.4 信號(hào)能量權(quán)重

電力系統(tǒng)低頻振蕩大多是多個(gè)模式的組合，其中阻尼比小且振蕩幅度大的主導(dǎo)模式特征信息是電力工作者最為關(guān)心的。當(dāng)信噪比足夠大時(shí)以信號(hào)能量作為定性分析指標(biāo)，從EEMD結(jié)果中經(jīng)過(guò)濾波去偽后提取真實(shí)模式中能量權(quán)重最大的一個(gè)IMF分量作為主導(dǎo)模式加以分析。

信號(hào)x的能量表達(dá)式為：

其中，t0、t1分別為仿真開(kāi)始和結(jié)束的時(shí)刻。

信號(hào)能量權(quán)重比定義為：

其中，Ef（i）為真實(shí)模式中第 i個(gè) IMF 分量的能量；m為真實(shí)模式IMF分量的總個(gè)數(shù)。

2 NExT

NExT是由JAMES等提出的一種適合于環(huán)境激勵(lì)的時(shí)域模式參數(shù)辨識(shí)方法，它的基本思想為：線性系統(tǒng)在白噪聲環(huán)境激勵(lì)下，結(jié)構(gòu)中2點(diǎn)之間響應(yīng)的互相關(guān)函數(shù)和脈沖響應(yīng)函數(shù)有相似的表達(dá)式，求得互相關(guān)函數(shù)之后就可以運(yùn)用時(shí)域模式辨識(shí)方法對(duì)其進(jìn)行辨識(shí)。

對(duì)于自由度為Q的線性系統(tǒng)，于系統(tǒng)k點(diǎn)處加激勵(lì) fk（t），則系統(tǒng) i點(diǎn)的響應(yīng) xik（t）為：

其中，φir為第i測(cè)點(diǎn)的第r階模式振型；akr為僅與激勵(lì)點(diǎn)k和模式階次r有關(guān)的常數(shù)項(xiàng)；λr為系統(tǒng)的第r階特征值。

根據(jù)振動(dòng)模式理論，當(dāng)系統(tǒng)的k點(diǎn)受到單位脈沖激勵(lì)時(shí)，則系統(tǒng)i點(diǎn)的脈沖響應(yīng)函數(shù)hik（t）為：

當(dāng)系統(tǒng)k點(diǎn)處受到激勵(lì)fk（t）時(shí)，系統(tǒng)i點(diǎn)和j點(diǎn)的響應(yīng) xik（t）和 xjk（t）的互相關(guān)函數(shù)為：

其中，E[·]表示求期望。

假定激勵(lì)f（t）是理想白噪聲，根據(jù)相關(guān)函數(shù)的定義，則有：

其中，δ為單位脈沖激勵(lì)；ak為僅與激勵(lì)點(diǎn)k有關(guān)的常數(shù)項(xiàng)。

將式（12）代入式（11）并積分，得：

對(duì)式（13）的積分部分進(jìn)行計(jì)算并化簡(jiǎn)，得：

將式（14）代入式（13），得：

對(duì)式（15）做進(jìn)一步的化簡(jiǎn)，經(jīng)整理得：

其中，bjr為僅與參考點(diǎn)j和模式階次r有關(guān)的常數(shù)項(xiàng)。

對(duì)比式（16）和式（10），可以發(fā)現(xiàn)兩者的數(shù)學(xué)表達(dá)式在形式上是完全一致的，因此互相關(guān)函數(shù)具有和系統(tǒng)的脈沖響應(yīng)函數(shù)同樣的性質(zhì)，可與其他模式辨識(shí)方法結(jié)合起來(lái)進(jìn)行環(huán)境激勵(lì)下的模式識(shí)別。

3 低頻振蕩主導(dǎo)模式識(shí)別

3.1 模式識(shí)別

3.1.1 Teager能量算子

Teager能量算子[16-18]具有健壯性和快速響應(yīng)能力，能迅速跟蹤信號(hào)的幅值和頻率。

設(shè)幅值和頻率均時(shí)變的工程信號(hào)為：

由文獻(xiàn)[17]可知信號(hào)xn的瞬時(shí)數(shù)字角頻率、實(shí)際頻率、幅值分別為：

EEMD得到的IMF分量頻率變化非常小，由于邊際效應(yīng)兩端數(shù)據(jù)存在誤差，計(jì)算過(guò)程中取合適的中間段數(shù)據(jù)求均值。

3.1.2 能量分析法

信號(hào)能量分析法[19]的原理是基于阻尼耗散能量思想，可將阻尼特性從用半周期能量描述的物理過(guò)程中提取出來(lái)。

振蕩信號(hào)x（t）的半周期能量定義為：

其中，ti為某一過(guò)零點(diǎn)時(shí)刻；Td為振蕩周期。

文獻(xiàn)[19]已給出信號(hào) x（t）的阻尼比 ξ為：

3.1.3 峰值法

由文獻(xiàn)[20]可知，已知頻率f和阻尼比ξ后，相位θ為：

其中，T（j）為信號(hào)的最大正峰值時(shí)間；S為最大正峰值個(gè)數(shù)；-π≤θ≤π。

將 j=1，2，…，S 代入式（24），可得：

由于-π≤θ≤π，可得：

由于g取整數(shù)，故可由式（26）確定g值，然后根據(jù)式（25）得到相位角。

3.2 基于EEMD-NExT的主導(dǎo)模式識(shí)別與預(yù)警

基于EEMD-NExT的低頻振蕩主導(dǎo)模式工況在線辨識(shí)與預(yù)警流程如圖2所示，具體步驟如下：

a.通過(guò)WAMS平臺(tái)獲取系統(tǒng)正常運(yùn)行下的2個(gè)信號(hào)并對(duì)其進(jìn)行預(yù)處理，然后分別作EEMD；

b.對(duì)EEMD得到的IMF分量由EEMD濾波器設(shè)置閾值獲取低頻振蕩范圍內(nèi)的分量；

c.利用互相關(guān)系數(shù)識(shí)別所得分量的真假，將虛假模式予以剔除，得到真實(shí)模式；

圖2 基于EEMD-NExT的低頻振蕩主導(dǎo)模式辨識(shí)與預(yù)警Fig.2 Identification and warning of low-frequency oscillation dominant mode based on EEMD-NExT

d.借助能量權(quán)重比對(duì)所得的真實(shí)模式分量進(jìn)行排序，選擇能量權(quán)重比最大的分量作為主導(dǎo)模式分量；

e.通過(guò)NExT對(duì)2個(gè)主導(dǎo)模式分量求互相關(guān)函數(shù)，將其作為主導(dǎo)模式辨識(shí)的信號(hào)；

f.對(duì)所得互相關(guān)函數(shù)分別通過(guò)Teager能量算子求時(shí)變幅值、時(shí)變頻率，通過(guò)能量分析求取阻尼比，時(shí)域峰值法求相位；

g.比較主導(dǎo)模式阻尼比與預(yù)警閾值大小，判斷系統(tǒng)是否發(fā)出告警。

4 算例仿真

4.1 數(shù)值信號(hào)算例

取復(fù)合數(shù)值信號(hào)

驗(yàn)證EEMD的抗模式混疊效果，以及依據(jù)相關(guān)系數(shù)甄別虛假模式的有效性和本文模式識(shí)別方法的可行性。

信號(hào)xx的EEMD結(jié)果如圖3所示。設(shè)定EEMD時(shí)空濾波器的閾值，此處設(shè)置帶通濾波頻率范圍為0.1～2 Hz，得到 imf7、imf8、imf9、imf10 這 4 個(gè) IMF 分量，求取它們與原信號(hào)xx的相關(guān)系數(shù)如表1所示，并根據(jù)本文模式識(shí)別方法辨識(shí)各分量的模式特征，結(jié)果如表2所示。

由表1可知，imf7和imf10的互相關(guān)系數(shù)均小于0.4；imf8和imf9的互相關(guān)系數(shù)分別為0.6469、0.9152，均比較大。由此判定imf7和imf10分量為偽分量。比較表2中imf7、imf8的模式信息和表3模式2的理論值可知，盡管兩者頻率很接近，但imf8的其他模式信息更接近理論值，由此進(jìn)一步驗(yàn)證了imf 7為虛假模式；表2中imf10的頻率為0.1638Hz，真實(shí)模式中并不存在該分量，可見(jiàn)根據(jù)互相關(guān)系數(shù)判斷imf10為虛假分量是有效的。表3為2個(gè)真實(shí)模式的理論值和本文方法的辨識(shí)結(jié)果，本文方法辨識(shí)誤差較小，是一種低頻振蕩模式辨識(shí)的有效方法。

4.2 EPRI-36節(jié)點(diǎn)算例

本文選取電科院的EPRI-36節(jié)點(diǎn)系統(tǒng)作為仿真算例，驗(yàn)證本文方法在多機(jī)系統(tǒng)研究中的有效性。

為了比較真實(shí)地模擬系統(tǒng)工況運(yùn)行的特點(diǎn)，本文選擇負(fù)荷隨機(jī)擾動(dòng)作為環(huán)境激勵(lì)，分別在負(fù)荷9、19和20上設(shè)置幅度為0.2 p.u.的功率隨機(jī)波動(dòng)。文獻(xiàn)[21]表明通過(guò)等值單機(jī)曲線可以很好地提取出主導(dǎo)振蕩信息，為此根據(jù)擴(kuò)展等面積準(zhǔn)則（EEAC）理論對(duì)系統(tǒng)進(jìn)行分群得，機(jī)組 1、2、3、4、5、7、8 為 S 群，機(jī)組6為A群，分別求取S群和A群的慣量中心δS、δA，并對(duì)它們添加λSNR=20 dB的高斯白噪聲以盡可能地再現(xiàn)工程實(shí)際中信號(hào)受噪聲干擾的特點(diǎn)。分別對(duì)加了高斯白噪聲的 δS和 δA進(jìn)行 EEMD，此處限于篇幅僅給出δS的EEMD效果圖，如圖4所示。

圖3 xx及EEMD效果圖Fig.3 xxand its EEMD results

表1 EEMD時(shí)空濾波后的IMF分量互相關(guān)系數(shù)Table 1 Cross-correlation coefficient of IMF components after EEMD spatiotemporal filtering

表2 EEMD時(shí)空濾波后的IMF分量模式信息Table 2 Mode information of IMF components after EEMD spatiotemporal filtering

表3 本文方法辨識(shí)結(jié)果與理論值比較Table 3 Comparison between results identified by proposed method and theoretical values

圖4 δS的 EEMD 效果圖Fig.4 δSand its EEMD results

區(qū)間振蕩模式較本地振蕩模式危害更大，其頻率范圍為0.1～1Hz，將該范圍留足裕度后設(shè)置為EEMD時(shí)空濾波器的帶通范圍，求取濾波后的各IMF分量的互相關(guān)系數(shù)見(jiàn)表4。對(duì)imf6、imf7求取信號(hào)能量后，并計(jì)算其能量權(quán)重比從而可得能量權(quán)重排序，結(jié)果如表4所示，可知imf6為主導(dǎo)模式。

取δS和δA中對(duì)應(yīng)的主導(dǎo)模式通過(guò)NExT求取互相關(guān)函數(shù)如圖5所示，再對(duì)其通過(guò)Teager能量算子求取主導(dǎo)模式頻率，通過(guò)能量分析求阻尼比。表5列出了本文方法、EMD-TEO-能量分析算法[18]、Prony 算法、小干擾分析主導(dǎo)模式辨識(shí)結(jié)果，其中EMD-TEO-能量分析算法和Prony算法的辨識(shí)結(jié)果均為對(duì)發(fā)電機(jī)7與發(fā)電機(jī)1的相對(duì)功角信號(hào)加噪20 dB后辨識(shí)所得，采樣頻率均為100 Hz。

表4 EEMD時(shí)空濾波后的IMF分量互相關(guān)系數(shù)Table 4 Cross-correlation coefficient of IMF components after EEMD spatiotemporal filtering

圖5 NExT法得到的互相關(guān)函數(shù)Fig.5 Cross-correlation function by NExT method

表5 4種方法對(duì)含噪信號(hào)辨識(shí)結(jié)果對(duì)比Table 5 Comparison of identification results among four methods for signal with noise

由表5可知，當(dāng)噪聲為20 dB時(shí)，4種方法均能辨識(shí)出EPRI-36節(jié)點(diǎn)系統(tǒng)的主導(dǎo)模式，本文方法和EMD-TEO-能量分析的辨識(shí)結(jié)果精度均高于傳統(tǒng)的Prony算法。本文方法的辨識(shí)結(jié)果更接近理論值，具有良好的抗噪性能，相比其他辨識(shí)方法無(wú)需人工激勵(lì)，是一種能有效在線辨識(shí)工況模式的方法。一般認(rèn)為，機(jī)電振蕩模式的阻尼比小于0.03時(shí)系統(tǒng)會(huì)發(fā)生低頻振蕩失穩(wěn)，故本文設(shè)置預(yù)警閾值為0.03，由辨識(shí)得到的主導(dǎo)模式阻尼比在閾值以下，此時(shí)系統(tǒng)發(fā)出告警須采取有效控制措施。

5 結(jié)論

a.EEMD與傅里葉變換、Prony分析等傳統(tǒng)辨識(shí)方法相比，更適合處理實(shí)際工程中的非線性、非平穩(wěn)工況信號(hào)，無(wú)需考慮定階問(wèn)題，同時(shí)抗噪能力強(qiáng)且能改善模式混疊現(xiàn)象。

b.通過(guò)EEMD時(shí)空濾波器對(duì)IMF分量濾波可以得到低頻振蕩工作范圍內(nèi)的信號(hào)，互相關(guān)系數(shù)能將虛假模式從真實(shí)模式中辨識(shí)出來(lái)并將其從中剔除，能量權(quán)重排序能找出主導(dǎo)模式分量。

c.本文采用了NExT，它是環(huán)境激勵(lì)模式識(shí)別的有效方法，避免了傳統(tǒng)模式辨識(shí)方法需人工激勵(lì)給電力系統(tǒng)安全帶來(lái)的危害，辨識(shí)出的模式參數(shù)更符合實(shí)際情況。

d.本文方法可以實(shí)時(shí)快速地跟蹤提取系統(tǒng)的主導(dǎo)模式信息，因此其在基于WAMS實(shí)測(cè)數(shù)據(jù)的低頻振蕩分析、在線監(jiān)測(cè)、預(yù)警及阻尼控制器設(shè)計(jì)等方面具有較好的實(shí)際應(yīng)用價(jià)值。