基于MEEMD與MP的電網(wǎng)低頻振蕩模態(tài)識(shí)別

張 程，劉佳靜，蔡思靜，林谷青，匡 宇

（1.福建工程學(xué)院電子電氣與物理學(xué)院，福建福州 350118；2.智能電網(wǎng)仿真分析與綜合控制福建省高校工程研究中心，福建福州 350118）

0 引言

電力系統(tǒng)低頻振蕩已經(jīng)成為互聯(lián)電網(wǎng)中的一大安全隱患，振蕩一旦發(fā)生便會(huì)持續(xù)較長一段時(shí)間才會(huì)消失，或繼續(xù)加劇誘發(fā)相繼事故，因此及時(shí)并準(zhǔn)確辨識(shí)出低頻振蕩信號(hào)模態(tài)參數(shù)具有重要意義。

目前用在電力系統(tǒng)低頻振蕩模式識(shí)別的方法有基于模型解析的特征根方法、傅里葉分析、旋轉(zhuǎn)不變性參數(shù)估計(jì)技術(shù)（Estimation of Signal Parameters Via Rotational Invariance Technique，ESPRIT）及其改進(jìn)算法、Prony 算法及Hilbert-Huang 變換等[1-8]。但是模型求取振蕩模式無法計(jì)及時(shí)變因素的影響；傅里葉分析基本只能分析穩(wěn)定信號(hào)，且完全喪失時(shí)域信息；希爾伯特黃變換（Hilbert-Huang Transform，HHT）雖能較好地反映信號(hào)的非平穩(wěn)性，但其在工程應(yīng)用中自身會(huì)存在虛假分量、能量泄露和端點(diǎn)效應(yīng)等問題；ESPRIT 算法在純信號(hào)的參數(shù)提取上較好，但在低信噪比下的性能較低且在模態(tài)參數(shù)估計(jì)時(shí)并沒有將復(fù)觀測(cè)數(shù)據(jù)及其共軛數(shù)據(jù)同時(shí)利用[9-10]；近年來，矩陣束（Matrix Pencil，MP）方法在此方向中被研究及應(yīng)用，文獻(xiàn)[11-12]采用MP 法來對(duì)振蕩及電路電流進(jìn)行在線辨識(shí)，但該方法受外在噪聲的干擾會(huì)出現(xiàn)較多虛假模態(tài)從而影響特征參數(shù)的提取。

在工程中一般會(huì)將量測(cè)數(shù)據(jù)先進(jìn)行預(yù)處理，較為常用的有自適應(yīng)濾波、卡爾曼濾波[13]和小波變換。但卡爾曼濾波需要預(yù)先確定系統(tǒng)模型，難以得到實(shí)際輸入信號(hào)特性；自適應(yīng)濾波雖算法簡單，但步長太大容易發(fā)散，步長太小雖是穩(wěn)定的但收斂速度慢[14]；小波去噪法雖計(jì)算量較小，但其存在閾值估計(jì)困難的限制。近些年來，有學(xué)者提出一種改進(jìn)的集成經(jīng)驗(yàn)?zāi)B(tài)分解法（Modified Ensemble Empirical Mode Decomposition，MEEMD），并受到了很大重視。它采用排列熵（Permutation Entropy，PE）[15]檢測(cè)信號(hào)的隨機(jī)性用于噪聲的消除，且可以有效改善經(jīng)驗(yàn)?zāi)B(tài)分解（Empirical Mode Decomposition，EMD）易出現(xiàn)模態(tài)混疊的現(xiàn)象。文獻(xiàn)[16-17]發(fā)現(xiàn)MEEMD能更清楚表達(dá)出鐵路故障信號(hào)的時(shí)頻分布。文獻(xiàn)[18]用MEEMD 實(shí)現(xiàn)人體生命體征的檢測(cè)。但是該新興方法目前還未在電力系統(tǒng)低頻振蕩方向有所應(yīng)用，本文將此法首先應(yīng)用在低頻振蕩模態(tài)辨識(shí)上。

本文提出MEEMD 與矩陣束相結(jié)合的方法并將其應(yīng)用于電網(wǎng)中低頻振蕩的模態(tài)辨識(shí)。該新方法首先對(duì)量測(cè)信號(hào)進(jìn)行MEEMD 分解獲得一系列本征模態(tài)函數(shù)（Intrinsic Mode Function，IMF），而后對(duì)除剩余殘量外的IMF 分量進(jìn)行重構(gòu)達(dá)到降噪目的，最后將重構(gòu)信號(hào)作為新的輸入信號(hào)進(jìn)行矩陣束分析提取低頻振蕩各模態(tài)特征。通過仿真驗(yàn)證本文提出的方法能夠有效抑制模態(tài)混淆并完整準(zhǔn)確辨識(shí)振蕩模態(tài)，相較于其他方法具有一定的抗噪性和優(yōu)越性。

1 MEEMD基本原理

對(duì)于非平穩(wěn)原始信號(hào)s(t)，在其中加入成對(duì)的振幅和標(biāo)準(zhǔn)差相等且方向相反的白噪聲ni(t) 和-ni(t)，得到的噪聲信號(hào)為：

式中：ai為添加的白噪聲的幅值；h為添加的白噪聲的對(duì)數(shù)。

根據(jù)大量實(shí)踐表明添加的幅值一般為原始信號(hào)標(biāo)準(zhǔn)差的10%～20%；h一般小于100 就可達(dá)到效果。

使用經(jīng)驗(yàn)?zāi)B(tài)分解分別對(duì)添加了高斯白噪聲的信號(hào)si+(t)和si-(t)進(jìn)行處理得到數(shù)個(gè)本征模態(tài)分量Iij+(t)和Iij-(t)，這樣信號(hào)可表示為：

式中：n為分解的IMF 個(gè)數(shù)；r（t）為余量。

將生成IMF 按式（3）進(jìn)行平均計(jì)算得到的結(jié)果記為Ij(t)，表示為：

將平均后的IMF 每個(gè)進(jìn)行排列熵的計(jì)算，結(jié)果記為(m)：

式中：m為嵌入維數(shù)；l為每個(gè)向量的序號(hào)；k為分解的總數(shù)；為每個(gè)分解量發(fā)生的概率。

當(dāng)=1/m!時(shí)，(m)≈ln(m!)，此時(shí)達(dá)到最大值，可以將排列熵(m)進(jìn)行標(biāo)準(zhǔn)化處理，即：

因此0 ≤≤1，由式（5）可知與時(shí)間序列的隨機(jī)程度成負(fù)相關(guān)，而噪聲信號(hào)的隨機(jī)性與排列熵值成正相關(guān)。信號(hào)的隨機(jī)性小，當(dāng)分解出的序列分量為平穩(wěn)信號(hào)時(shí)，排列熵值就較小[19]，因此需要設(shè)置閾值θ來區(qū)分信號(hào)的熵值用于干擾信號(hào)的剔除，經(jīng)過多次實(shí)驗(yàn)，發(fā)現(xiàn)0.5～0.6 比較合適，所以本文設(shè)置閾值θ=0.6。當(dāng)分解出的序列分量的熵值大于閾值θ時(shí)，將其視為異常信號(hào)，這樣依次判斷，最后將異常信號(hào)累加。

篩選完畢后從原始信號(hào)s(t)將所有累加異常信號(hào)剔除后再進(jìn)行EMD 分解。所以最后MEEMD 分解表達(dá)式為：

式中：r(t)為最終得出的剩余殘量，在后續(xù)的重構(gòu)分析中需將此分量剔除。

2 矩陣束分析

設(shè)信號(hào)形式為：

式中：M為多個(gè)模態(tài)函數(shù)的階數(shù)；y(t)為響應(yīng)信號(hào)；αi為衰減因子；fi為振蕩頻率；ε(t)為白噪聲信號(hào)；y0(t)為實(shí)際信號(hào)。

其中，

式中：Ai為幅值；θi為初始的相位。

對(duì)于確定最大階數(shù)M可以通過采樣序列構(gòu)造HANKEL 矩陣，如下所示：

式中：N為采樣數(shù)的最大值；L為MP 參數(shù)，通常N/4 ≤L≤N/3。

對(duì)Y矩陣進(jìn)行奇異值分解[20]，即：

式中：U為(N-L)×(N-L) 階的正交矩陣；V為(L+1)×(L+1)階的正交矩陣；S為(N-L)×(L+1)階的對(duì)角陣。

可以通過設(shè)置適當(dāng)?shù)拈撝祦硖崛∏癕個(gè)奇異值的最大階數(shù)[21]，從而形成1 個(gè)新的矩陣S1：

取矩陣V的前2M個(gè)奇異向量構(gòu)成V1，構(gòu)造如下2 個(gè)矩陣：

式中：V11和V12分別為刪除V1的最后一行和第一行得到的新矩陣。

定義zi為Y2-λY1（λ為約束值）的廣義特征值，然后使用式（13）便可得每個(gè)振蕩模式的頻率fi和衰減因子αi，即：

式中：Re()· 與Im()· 分別表示zi的實(shí)部與虛部；Δt為采樣間隔時(shí)間。

3 評(píng)價(jià)指標(biāo)

3.1 信噪比

為評(píng)價(jià)所提方法的去噪效果，以信噪比（Signal to Noise Ratio，SNR）為測(cè)量指標(biāo)[22]，計(jì)算公式為：

式中：Ps為信號(hào)功率；Pn為噪聲功率。

信噪比越高，信號(hào)越純凈，混合噪聲成分越少。

3.2 相關(guān)系數(shù)

為了評(píng)價(jià)分量序列重構(gòu)后信號(hào)與原始信號(hào)的匹配度，這里采用相關(guān)系數(shù)Coef來作為指標(biāo)[23]，其具體的計(jì)算公式為：

式中：Xi為原始信號(hào)序列，i=1，2，...，N；Yi為重構(gòu)信號(hào)序列，i=1，2，...，N；N為采樣點(diǎn)數(shù)。

Coef的區(qū)間為[-1，1]，其絕對(duì)值與相關(guān)性為正比關(guān)系，|Coef|越逼近于1，則相關(guān)性越強(qiáng)，反之則越弱。

4 仿真算例

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

為驗(yàn)證本文方法的有效性，構(gòu)造1 個(gè)低頻振蕩信號(hào)。假設(shè)：

該信號(hào)包含了2 個(gè)頻率分別為0.7 Hz 和1 Hz的低頻振蕩模式稱為模態(tài)1、模態(tài)2。其衰減因子分別為-0.1 和-0.5，v(t)為信噪比為9.5 dB 的高斯白噪聲。由于這2 個(gè)模態(tài)頻率非常相近，只相差了0.3 Hz，且里面混有噪聲干擾，這必然會(huì)引起模態(tài)混疊現(xiàn)象。通過此來驗(yàn)證本文所提方法在低頻振蕩并伴有模態(tài)混疊情況下的辨識(shí)準(zhǔn)確度。采樣頻率為25 Hz，采樣點(diǎn)數(shù)為500，原始信號(hào)為圖1 所示，其中，A為幅值，無量綱。

圖1 原始信號(hào)Fig.1 Original signal

圖2 所示的一系列IMF 分量是低頻振蕩信號(hào)分別經(jīng)EMD 和MEEMD 分解后所得。

圖2 低頻振蕩信號(hào)被分解后的IMF分量Fig.2 IMF component after the oscillation signal is decomposed

由圖2 可知，在低信噪比的情況下MEEMD 只產(chǎn)生了4 個(gè)本征模態(tài)函數(shù)和1 個(gè)剩余殘量，而EMD分解出了6 個(gè)本征模態(tài)函數(shù)和1 個(gè)殘余量，且EMD分解后的第1 個(gè)IMF 分量和第2 個(gè)IMF 分量由于受到噪聲干擾嚴(yán)重而頻率混淆波動(dòng)，受到污染，相當(dāng)于無用虛假分量，但MEEMD 分解并無此現(xiàn)象?？梢奙EEMD 能夠抑制噪聲干擾和模式混疊，減少虛假分量。

圖3 所示的2 組信號(hào)分別為含噪信號(hào)經(jīng)EMD與MEEMD 分解后再過濾掉剩余殘量后的重構(gòu)信號(hào)，其各方面指標(biāo)見表1。

圖3 2種方法分解后的重構(gòu)信號(hào)Fig.3 Reconstructed signals decomposed by two methods

表1 重構(gòu)信號(hào)各指標(biāo)對(duì)比Table 1 Index comparison of reconstructed signal

通過對(duì)比可以發(fā)現(xiàn)信號(hào)經(jīng)MEEMD 處理后可以將信噪比為原來的9.5 dB 提升到20.43 dB，且均方誤差也小于EMD，與原始信號(hào)的相關(guān)系數(shù)高達(dá)0.993 5，再次證明本文使用MEEMD 作為低頻振蕩信號(hào)預(yù)處理的方法在抑制模態(tài)混淆和抗噪性上具有一定的優(yōu)勢(shì)。

為了說明本文方法的可行有效性，將低頻振蕩信號(hào)分別用本文方法、傳統(tǒng)HHT 法和MP 法進(jìn)行處理，辨識(shí)結(jié)果如表2 所示。

表2 振蕩信號(hào)經(jīng)多種方法辨識(shí)結(jié)果Table 2 Oscillation signal identified by various methods

由表2 參數(shù)對(duì)比可知，在噪聲干擾且頻率易混疊的情況下，本文方法都能夠準(zhǔn)確提取出各個(gè)模態(tài)的特征參數(shù)且誤差較?。欢溆? 種方法由于噪聲干擾誤差較大，且HHT 因?yàn)榛殳B現(xiàn)象，甚至辨識(shí)出了幾個(gè)高頻虛假分量。所以本文方法在低信噪下不僅能夠準(zhǔn)確提取低頻振蕩特征參數(shù)，且誤差較小，再次說明其抗噪性和精確性。

4.2 EPRI-36節(jié)點(diǎn)系統(tǒng)

本文采用PSASP 中EPRI-36 節(jié)點(diǎn)系統(tǒng)[24-25]進(jìn)行分析，考慮如下擾動(dòng)：在1 s 時(shí)BUS20 至BUS22 之間聯(lián)絡(luò)線20%處發(fā)生持續(xù)時(shí)間為0.1 s 的三相短路故障，采樣步長為0.01 s，總時(shí)長為15 s。采集發(fā)電機(jī)G3 相對(duì)功角信號(hào)（以G1 為參考發(fā)電機(jī)）作為分析振蕩信號(hào)，并在該信號(hào)中人工添加信噪比為9.5 dB的高斯白噪聲得到的含噪振蕩信號(hào)，如圖4 所示。

圖4 G 3功角振蕩信號(hào)Fig.4 G3 power angle oscillation signal

將含噪信號(hào)分別通過EMD 和MEEMD 分解得到的一系列IMF 分量如圖5 所示。

圖5 振蕩信號(hào)分別被分解后的IMF分量Fig.5 IMF component after oscillation signal is decomposed

從圖5（b）中可知，EMD 分解出來的第1 個(gè)IMF分量由于受到噪聲干擾而頻率嚴(yán)重混淆，影響后期模態(tài)辨識(shí)精度；但在相同噪聲干擾下MEEMD 并無此況且分解的IMF 分量個(gè)數(shù)較少，可見此法能有效降低干擾成分，抑制模態(tài)混疊并減少虛假成分。

圖6 給出了振蕩信號(hào)經(jīng)分解后的重構(gòu)信號(hào)，各方面的指標(biāo)見表3。

圖6 重構(gòu)信號(hào)Fig.6 Reconstructed signals

表3 重構(gòu)信號(hào)各指標(biāo)對(duì)比Table 3 Index comparison of reconstructed signal

由表3 可知，經(jīng)MEEMD 處理后重構(gòu)信號(hào)可將信噪比由原來的9.5 dB 提升到20.263 5 dB，可見其降噪效果的優(yōu)勢(shì)，且在與原信號(hào)的相關(guān)性和均方差指標(biāo)上均優(yōu)于另一種。

分別采用MP，HHT 和本文方法對(duì)含噪功角振蕩信號(hào)進(jìn)行模態(tài)辨識(shí)，結(jié)果見表4。

表4 多種方法識(shí)別結(jié)果Table 5 Recognition result by multiple methods

由表4 參數(shù)對(duì)比可知，在噪聲干擾下，本文方法能夠辨識(shí)出2 種頻率，并與小干擾穩(wěn)定計(jì)算結(jié)果（0.777 5 Hz 和0.980 2 Hz）相吻合；但是在相同情況下，MP 法在提取參數(shù)特征的精度上跟本文方法比還是稍有不足的，HHT 甚至無法辨識(shí)出另一個(gè)模態(tài)?？梢姳疚奶岢龅姆椒ㄔ诙鄼C(jī)多點(diǎn)的互聯(lián)大電網(wǎng)中的低頻振蕩模態(tài)辨識(shí)也是準(zhǔn)確有效的。

5 結(jié)論

本文提出MEEMD-MP 新方法用于低頻振蕩模態(tài)辨識(shí)，仿真得到以下結(jié)論：

1）考慮到實(shí)際電網(wǎng)低頻振蕩信號(hào)的特征，利用MEEMD 作為數(shù)據(jù)預(yù)處理手段，同時(shí)根據(jù)噪聲信號(hào)和源信號(hào)隨機(jī)性的不同引入排列熵用于篩選和分離，既能在保留原始信號(hào)特征的情況下抑制模態(tài)混疊，降低噪聲干擾，提高信噪比。

2）將分解后并重構(gòu)的信號(hào)作為新的輸入信號(hào)通過矩陣束處理，改善了矩陣束方法因有噪聲干擾而無法獲取理想的低頻振蕩模態(tài)信息的不足。

3）仿真分析表明本文提出的方法具有一定的抗干擾性和精確性，是一種行之有效的電力系統(tǒng)低頻振蕩模式識(shí)別方法。