閆紅艷,Hwang Jin Kwon,高艷豐
(1. 河北工程大學水利水電學院,河北省 邯鄲市 056038;2. 韓國又石大學能源工程系,韓國鎮(zhèn) 川365-803)
隨著電網(wǎng)遠距離、大容量輸電的實施,高放大倍數(shù)勵磁裝置的使用以及大規(guī)模新能源的接入,低頻振蕩時有發(fā)生,頻率范圍一般為0.1~2.5 Hz[1],嚴重危害到系統(tǒng)正常運行。低頻振蕩模態(tài)識別是分析低頻振蕩產(chǎn)生原因的基礎(chǔ),也是實現(xiàn)電力系統(tǒng)實時高效控制和風險預警的關(guān)鍵之一[2-5]。
基于相量測量單元(phasor measurement unit,PMU) 的電力系統(tǒng)廣域測量系統(tǒng)(wide area monitoring systems,WAMS)能夠?qū)崟r地記錄電力系統(tǒng)的動態(tài)行為,按照量測信號的辨識方法可分為基于大擾動后自由振蕩響應信號的方法和基于環(huán)境激勵下隨機響應信號的方法[6-10]。大擾動數(shù)據(jù)包含強低頻振蕩分量,能夠準確識別區(qū)域間模式,但這類事件很少發(fā)生,很難進行連續(xù)的模態(tài)識別。相反,負荷隨機波動引起的小擾動激勵(又稱類噪聲數(shù)據(jù))時刻存在、易于采集、數(shù)據(jù)豐富,可及時準確地反映當前系統(tǒng)的運行特性[11]。Prony 算法及其改進算法在分析大擾動激勵下的響應信號應用較廣泛,但Prony 算法在處理噪聲信號和具有時變特性非平穩(wěn)信號的能力不理想[12-13]。基于環(huán)境激勵下隨機響應信號的方法大多在信號的自回歸(auto regressive,AR)模型[10]或自回歸滑動平均(auto regressive moving average,ARMA)模型上發(fā)展起來的[14-16]。文獻[17-18]采用最小二乘法求取電力系統(tǒng)隨機響應信號AR 模型的最優(yōu)參數(shù),進而獲得低頻振蕩的頻率和阻尼比信息,但是該算法存在數(shù)值不穩(wěn)定的問題。文獻[19]提出了用Yule-Walker 方程解AR 模型的功率譜分析方法;文獻[20]提出了一種基于數(shù)學形態(tài)學自回歸移動平均(mathematical morphology,MM-ARMA)算法的辨識方法。文獻[21]釆用Bootstrap 法確定YW 法辨識結(jié)果的置信區(qū)間。改進擴展的Yule-Walker(modified extended Yule-Walker,MEYW)方法[22]廣泛用于基于環(huán)境數(shù)據(jù)自相關(guān)函數(shù)的AR 模型中,但AR 模型的階數(shù)較難確定,存在2 個頻率相近的振蕩模式無法區(qū)分辨識的現(xiàn)象。
隨機子空間方法(stochastic subspace identification,SSI)從原理上既適用于自由振蕩信號分析,也適用于由環(huán)境激勵引起的系統(tǒng)隨機響應信號分析,抗噪性強,模態(tài)識別結(jié)果較準確,但計算速度較慢,很難反映信號的時變特性[7,23]。低頻振蕩產(chǎn)生的直流分量會降低模態(tài)識別的準確性,上述方法在使用過程中數(shù)據(jù)均需進行去除直流預處理。希爾伯特-黃變換(Hilbert-Huang transform,HHT)算法是分析非線性、非平穩(wěn)功率振蕩的常用工具,雖然能分解出直流分量,但存在端點效應現(xiàn)象[24]和本征模函數(shù)篩選效果不理想等固有缺點,辨識結(jié)果誤差率較大。文獻[25]提出集合經(jīng)驗模態(tài)分解(ensemble empirical mode decomposition,EEMD)算法,但該算法存在計算復雜、耗時長等問題。
變分模態(tài)分解(variational mode decomposition,VMD)是一種新的非遞歸的信號分解方法,除了能濾除直流外,還能克服經(jīng)驗模態(tài)分解(empirical mode decomposition,EMD)和EEMD 方法的缺陷,噪聲魯棒性較好。本文采用基于模態(tài)相關(guān)系數(shù)的VMD 進行類噪聲數(shù)據(jù)處理,直接去除中心頻率為0 的直流部分并分離出低頻振蕩信號;低頻振蕩信號模態(tài)參數(shù)辨識利用信號頻率差序列的自相關(guān)函數(shù)表示為與低頻振蕩函數(shù)形式一樣的指數(shù)衰減正弦函數(shù)的線性組合,該函數(shù)的固有頻率和阻尼比等于振蕩模態(tài)的固有頻率和阻尼比,用不同時間頻率數(shù)據(jù)之間的差序列減少測量噪聲的影響。通過在頻域中峰值頻率附近的離散傅里葉變換(discrete Fourier transform,DFT)進行曲線擬合來估計其拉普拉斯變換系數(shù),然后計算得出模態(tài)分量的參數(shù)。通過與MEYW方法的識別結(jié)果相比較,驗證了所提識別方法的準確性。最后,采用基于ARMA 模型產(chǎn)生PMU 頻率數(shù)據(jù)和某實際系統(tǒng)的實測頻率數(shù)據(jù)驗證方法的可行性和有效性。
VMD 是Dragomiretskiy 提出的非遞歸信號分解方法,實質(zhì)是變分問題,根據(jù)預設(shè)模態(tài)分量個數(shù)對信號進行分解[26]。該方法通過一個自適應維納濾波器組將原始信號f(x)分解為K個中心頻率為ωk的模態(tài)函數(shù)uk,其中K為預設(shè)模態(tài)分量個數(shù)。VMD算法在抗噪聲和非平穩(wěn)信號處理方面具有較好的性能和較高的運算效率,可以分解出直流分量。
為了得到具有一定帶寬頻率的K個模態(tài)分量,通常對每個模態(tài)函數(shù)uk進行Hilbert 變換得到邊際譜:

式中:δ(t)為狄利克雷函數(shù);j2=-1;*為卷積符號。預估各模態(tài)解析信號中心頻率,將每個模態(tài)的頻譜調(diào)制到相應的基頻帶:

式中:{uk}={u1,…,uK}為分解得到的K個模態(tài)分量;{ωk}={ω1,…,ωK}為各分量的頻率中心。
受約束的變分問題求解,可引入二次懲罰因子α和拉格朗日乘法算子λ(t),得到增廣拉格朗日公式:

利用乘法算子交替方向法求取式(6)變分問題,通過交替更新un+1k、ωn+1k和λn+1尋求增廣拉格朗日表達式的鞍點。式(5)為約束變分模型的最優(yōu)解,從而將f分解為K個窄帶IMF 分量。VMD 算法的具體過程如下:
1)初始化{u1k},{ω1k},λ^1,n=0;
2)n=n+1,執(zhí)行迭代循環(huán);
3)使k=k+1,按照式(5)與式(6)更新u^n+1k與ωn+1k,直至k=K;
4)按照式(7)更新λ^n+1;

VMD 本身不具有自適應性,模態(tài)分解數(shù)K值的設(shè)置是信號VMD 分解的關(guān)鍵環(huán)節(jié),對分解效果影響較大。K設(shè)定值大于待分解信號所含的固有模態(tài)(intrinsic mode function,IMF)分量個數(shù),則會在最終結(jié)果中引入虛假模態(tài)分量,影響對原始低頻振蕩信號的分析;相反,K設(shè)定值過小,則將導致信號分解不完全,即振蕩信號中含有的重要模態(tài)沒有被完全分解出來。
本文利用模態(tài)分量的相關(guān)系數(shù)來確定VMD分解個數(shù)的方法。首先取K值為2,計算各個分量之間的相關(guān)系數(shù),判斷各分量之間是否存在頻率混疊現(xiàn)象,自適應確定模態(tài)個數(shù)。分量x1(n)和分量x2(n)的相關(guān)系數(shù)定義為

基于模態(tài)相關(guān)系數(shù)的VMD 算法具體步驟如下:
1)初始化K= 2,并用VMD 算法處理原始信號;
2)計算各個模態(tài)分量之間的相關(guān)系數(shù),提取其中最大的相關(guān)系數(shù);
3)K→K+1,并根據(jù)步驟2)更新模態(tài)之間的最大相關(guān)系數(shù);
4)重復步驟3),直到最大相關(guān)系數(shù)超過閾值。經(jīng)過大量實驗分析,最大相關(guān)系數(shù)的閾值選取0.1較為合適[27]。
根據(jù)文獻[26],當懲罰參數(shù)α=2 000時,它可以滿足大多數(shù)工作條件的實際需求。試驗表明,改變α時影響很小,因此取默認值,當存在直流分量時,DC=1,這意味著可以分解出中心頻率為0 的分量。本文首要任務是要濾除直流分量,所以DC取1,其余參數(shù)取默認值。
負荷的隨機變化即使很小也會影響到供需之間的不平衡,通常將這種干擾看作高斯白噪聲[28],在電力系統(tǒng)類噪聲數(shù)據(jù)中用直流分量d(t)和低頻振蕩模式s(t)表示。同時,考慮到WAMS 數(shù)據(jù)中含有各種測量噪聲w(t),因此,電力系統(tǒng)正常運行時,PMU數(shù)據(jù)的信號模型可表示為

式中:J為系統(tǒng)的慣性常數(shù);f0為標稱頻率;Ks為系統(tǒng)的功率/頻率特性;u(t)為系統(tǒng)負載功率缺額的標幺值。d˙(t)是頻率波動引起的直流分量,即每日負荷波動時的功率缺額;u(t)是階躍函數(shù)的總和,是負荷切換時間的函數(shù)[30]。日常負荷隨機變化被建模為環(huán)境頻率數(shù)據(jù)中的白噪聲,所以d(t)以白噪聲v(t)為輸入時表示為

式中v(t)=-u˙(t),隨著v(t)的強度變強,由于電力系統(tǒng)中的調(diào)速器或頻率控制,隨著時間的推移,v(t)的強度可能會逐漸不同于其實際值。
擾動v(t)引起電力系統(tǒng)頻率中的低頻振蕩s(t)。振蕩模式被表示為二階微分方程[31]。K個振蕩信號可表示為

式中:w[n]是測量噪聲;d[n]=d(nTs);s[n]=s(nTs),Ts是PMU 的采樣周期。測量噪聲w[n]可以看成白色高斯信號,其均值和方差分別是0和σ2。y[n]的信噪比(signal noise ratio,SNR)定義為s[n]與w[n]的功率比。在模擬產(chǎn)生數(shù)據(jù)和模態(tài)識別時,可將式(11)和式(13)的連續(xù)時間方程轉(zhuǎn)換為離散時間方程,用脈沖響應不變法[32]進行轉(zhuǎn)換,表示成ARMA模型為:

將原始頻率信號y[n]進行VMD分解,分離出中心頻率為零的直流分量d[n],同時得到K-1個有限帶寬的固有本征模分量,根據(jù)其頻率值可直接提取出含噪聲的低頻振蕩分量sw[n],再對該模態(tài)分量進行參數(shù)辨識即可得系統(tǒng)的振蕩參數(shù)。含噪的低頻振蕩信號sw[n]的ARMA模型可以表示為

式中ck和ek分別為自回歸部分和滑動平均部分模型參數(shù)。
由信號的自相關(guān)特性可知,自相關(guān)分析能有效消除信號中的噪聲,且能保留原信號函數(shù)的頻率特征[33-34],隨著時間的延長,噪聲信號自相關(guān)函數(shù)值將很快衰減至0,因此,可以采用自相關(guān)函數(shù)代替原函數(shù)進行模態(tài)辨識。
故sw[n]的自相關(guān)函數(shù)可以表示為

rw[m]是與式(17)中的AR模型具有相同特征多項式的線性預測模型[34]。因此,rw[m]的連續(xù)信號也可以表示為K個指數(shù)衰減正弦分量的線性組合,表達式為

式中Re{·}代表復數(shù)的實數(shù)部分。
公式(19)自相關(guān)函數(shù)rw[m]進行DFT的結(jié)果為

為了減少曲線擬合中的噪聲影響和縮短計算時間,希望在fk附近使用Rw[m],fk是低頻振蕩模態(tài)功率最集中的地方。最好在以下頻率范圍上進行曲線擬合:|f^k-f|≤fc,k=1,2,…,K。
本文中將fc設(shè)為0.05 Hz,其寬度足以包含每個峰值的頻率范圍。每一個峰值曲線擬合采樣的最大整數(shù)小于2fc/△f。角頻率采樣點表示為ωk,m和ωk,(m+1)=ωk,m+2π△f。m=1,2,…,M是第K個模態(tài)峰值的曲線擬合。
信號VMD 分解后,取低頻振蕩頻率范圍內(nèi)IMF分量,并對此分量做自相關(guān)計算,按式(24)中R(s)的系數(shù)ak和bk來辨識模式信號各參數(shù)的值。整個過程通過將R(s)擬合到頻域中的R[m]中來實現(xiàn)。R(s)在頻域中可以表示為


式中?表示偽逆,X的解即式(24)的系數(shù)ak和bk。則第k個模態(tài)參數(shù)f^k,ζ^k,A^k和φ^k即可通過式(24)和式(25)計算得出。再由式(20)可得r[n]的擬合值r^[n]:

信號擬合曲線與原信號越接近,辨識結(jié)果精確度也越高。擬合精度采用信噪比SNRe為指標,單位為dB。

式中:r[n]為測量信號;r^[n]為曲線擬合重構(gòu)信號;SNRe的結(jié)果越大,表示擬合信號與原始信號擬合的效果越好。
本文采用VMD分解濾除原始數(shù)據(jù)中的直流分量,并提取出系統(tǒng)低頻振蕩信號,然后做自相關(guān)計算消除噪聲,低頻振蕩參數(shù)由自相關(guān)函數(shù)的拉普拉斯變換計算得到,而拉普拉斯變換式(22)的系數(shù)通過信號自相關(guān)函數(shù)DFT的曲線擬合來估計,最后由ARMA模型產(chǎn)生的類似真實PMU的頻率數(shù)據(jù)和實測PMU頻率數(shù)據(jù)驗證方法的有效性和準確性。
由式(9)可知,PMU數(shù)據(jù)由3部分組成,為了驗證所提方法的可行性和有效性,利用式(15)和式(16)的ARMA模型模擬產(chǎn)生PMU頻率數(shù)據(jù)進行驗證。由ARMA模型可通過程序模擬產(chǎn)生離散的PMU數(shù)據(jù),直流分量[27]取值如表1所示,2個低頻振蕩分量取值如表2 所示。模擬采樣時間間隔為0.02 s,加入不同信噪比的高斯白噪聲來模擬測量噪聲。

表1 模擬直流分量的參數(shù)Tab.1 Parameters of the analog DC component

表2 低頻振蕩分量的參數(shù)Tab.2 Parameters of low-frequency oscillation components
模擬數(shù)據(jù)測量噪聲信噪比取30 dB 高斯白噪聲時,對模擬數(shù)據(jù)先經(jīng)過低通濾波處理,低通濾波是基于高速采樣頻率50 Hz設(shè)計的2階巴特沃茲低通濾波器,截止頻率為5 Hz,濾除高頻分量可減少VMD 分解個數(shù),增加運行速度。采用相關(guān)系數(shù)的VMD算法確定分量個數(shù)K,分量之間的最大相關(guān)系數(shù)如表3所示。由表3可知,當K取4、5時,信號經(jīng)VMD算法分解之后分量之間的最大相關(guān)系數(shù)均大于閾值0.1;而當k取2、3時,最大相關(guān)系數(shù)皆小于閾值0.1。故取K值為3,VMD分解結(jié)果如圖1所示。

表3 不同K值的最大相關(guān)系數(shù)Tab.3 Maximum correlation coefficients of K different values
由圖1可知,信號經(jīng)VMD分解后的各個IMF呈現(xiàn)比較規(guī)范,彼此間沒有模態(tài)混疊現(xiàn)象,各個頻段分離效果較好。其中IMF0是中心頻率為0的直流分量,其余的是不同頻率范圍的主導振蕩頻率,為后續(xù)準確辨識出低頻振蕩特征參數(shù)提取提供了理想的模態(tài)分量。

圖1 VMD模態(tài)分量及頻譜Fig.1 The spectrum and VMD modal component
與EEMD算法相比,VMD算法具有較好的優(yōu)越性,EEMD 分解結(jié)果如圖2 所示。從圖中可以看出,經(jīng)過EEMD分解后得到13個IMF,右側(cè)為對應頻譜。分解得到的模態(tài)個數(shù)遠遠多于原始信號含有的振蕩分量個數(shù),且耗時很長,同時出現(xiàn)了模態(tài)混疊,無法準確反映原始信號的低頻振蕩分量,嚴重影響參數(shù)提取的準確度。
以圖1 中IMF1 和IMF2 信號為例進行分析,對IMF1 分別利用頻差序列DFT 的曲線擬合法和MEYW法進行振蕩分量模態(tài)辨識,辨識算法采樣頻率為50 Hz,自相關(guān)函數(shù)的有效持續(xù)時間設(shè)置為20 s,數(shù)據(jù)窗長為5 min,相鄰數(shù)據(jù)窗間隔為1 min,數(shù)據(jù)總長為15 min。得到3~8 min 的頻率偏差波形和20 s頻率偏差的自相關(guān)擬合曲線如圖3所示,擬合DFT 幅值和角度如圖4 所示。可以看出擬合曲線非常接近辨識信號的曲線,辨識結(jié)果準確度高。圖5給出IMF1分量的頻率和阻尼比的辨識結(jié)果,從圖中可以看出2 個方法的頻率偏差比較小,基于MEYW法的阻尼率值偏小,且受噪聲影響某個時間段誤差會比較大,而基于DFT的曲線擬合法變化比較平穩(wěn),抗噪性好。將各滑動窗口內(nèi)辨識得到的模態(tài)參數(shù)取平均值得到辨識結(jié)果。采用同樣的方法,將測量噪聲設(shè)成10 dB 高斯白噪聲時,K值為4,并對低頻振蕩頻率范圍內(nèi)的分量進行模態(tài)辨識,模擬數(shù)據(jù)辨識結(jié)果如表4所示。

圖3 模擬數(shù)據(jù)IMF1分量的頻率偏差波形及自相關(guān)函數(shù)Fig.3 Frequency deviation waveform and autocorrelation function of IMF1

圖4 模擬數(shù)據(jù)IMF1分量的DFT幅值和角度Fig.4 DFT amplitude and angle of IMF1

圖5 模擬數(shù)據(jù)IMF1頻率和阻尼率的辨識結(jié)果Fig.5 Frequencies and damping rates identification results of IMF1

表4 模擬數(shù)據(jù)的辨識結(jié)果Tab.4 Identification results of simulated data
由表4可以看出,考慮量測噪聲影響時,本文提出的辨識法能給出較為準確的頻率值,而阻尼比和振蕩幅值受系統(tǒng)運行方式的影響較大,所以信噪比越小,系統(tǒng)阻尼比和振蕩幅值波動程度相對越大,但本文提出的方法波動相對較小,說明所提方法準確度高和抗噪性好。因此,文中采用的基于DFT 曲線擬合的辨識法比目前應用廣泛的在線辨識MEYW法更準確,擬合精度SNRe如圖6所示。
對比圖6 中2 個分量的擬合精度SNRe可以看出,弱阻尼辨識結(jié)果比高阻尼辨識結(jié)果更準確,阻尼比越小,振蕩平息時間也越長,一段時間后該模式信號分量相對于噪聲仍占據(jù)主導地位,同時也可看出環(huán)境噪聲越弱,信噪比越高,辨識精度越高。

圖6 模擬數(shù)據(jù)擬合精度圖Fig.6 Fitting precision diagram
該方法的計算時間與VMD分解個數(shù)、采樣頻率、數(shù)序列長度、數(shù)據(jù)窗滑動步長以及時間窗長度直接相關(guān),VMD計算時間占比較大,但相對于EEMD 計算時間能減少很多,同時低通濾波后可減少分解個數(shù),提高運行速度。由于MEYW法中AR模型的階數(shù)比DFT擬合法的階數(shù)高,相比較,本文辨識方法運行時間短。以MATLAB 2018 版進行編程,完成上述PMU模擬數(shù)據(jù)分析,SNR為30 dB 時低通濾波后VMD 分解個數(shù)為3,計算時間為6.372 7 s,基于DFT 曲線擬合辨識運行時間為0.927 8 s,耗時7.518 9 s,占空比為0.835 4%,完全滿足在線應用要求。
以某實際系統(tǒng)PMU錄波數(shù)據(jù)為例,分析系統(tǒng)低頻振蕩特征。實際系統(tǒng)采樣頻率為0.033 3 s,選取時長為15 min 的部分量測信號,實測數(shù)據(jù)減去額定頻率60 Hz,得到實際頻率的波動如圖7所示。

圖7 某電網(wǎng)實測頻率波動Fig.7 The measured frequency fluctuation of a power grid
VMD分解時經(jīng)巴特沃斯低通濾波處理可以減少VMD分解個數(shù),提高運行速度。本案例中,模態(tài)個數(shù)為5時最大相關(guān)系數(shù)大于閾值0.1,故取模態(tài)個數(shù)為4。VMD分解得到各分量及其頻譜如圖8所示,得到第一個IMF 分量為中心頻率為0 的直流,同時得到3個不同頻率段的信號,取危害較大的區(qū)域間低頻振蕩的信號IMF1 和IMF3,用本文方法進行模態(tài)辨識,辨識時采樣頻率為30 Hz,自相關(guān)函數(shù)的有效持續(xù)時間為20 s,數(shù)據(jù)窗長設(shè)為5 min,相鄰數(shù)據(jù)窗間隔為1 min。此時得到辨識結(jié)果如表5 所示,IMF1分量的幅值和角度如圖9所示,實測數(shù)據(jù)2個區(qū)間低頻振蕩模態(tài)分量的頻率和阻尼比的辨識結(jié)果如圖10所示,擬合精度如圖11所示。通過對比結(jié)果可知,本文所提方法準確度高、抗噪性強,MEYW方法受隨機測量噪聲影響較大,同時通過對比擬合精度也驗證了阻尼率越小,擬合精度越高。

表5 實測數(shù)據(jù)辨識結(jié)果Tab.5 Identification results of measured data

圖8 某電網(wǎng)實測頻率數(shù)據(jù)的模態(tài)分量及頻譜Fig.8 Modal components and spectrum of a power network measured frequency data

圖9 實測數(shù)據(jù)IMF1分量的DFT幅值和相角Fig.9 IMF1 DFT amplitude and angle of measured data

圖10 實測數(shù)據(jù)2個模態(tài)的辨識結(jié)果Fig.10 Two modes identification results of measured data

圖11 實測數(shù)據(jù)擬合精度Fig.11 Fitting accuracy of measured data
在實測數(shù)據(jù)辨識過程中,采樣頻率較低,為30 Hz,同樣在MATLAB 2018 版本上運行程序。低通濾波后VMD 分解個數(shù)為4,運行時間為8.404 8 s,總耗時為9.551 s,占空比為1.061%,提高了機電小干擾穩(wěn)定評估的實時性。
以上結(jié)果表明,由于VMD 算法和基于DFT的曲線擬合法都具有較好的噪聲魯棒性,參數(shù)辨識結(jié)果精度較高,辨識曲線平滑,符合實際的波動情況,并且計算速度滿足在線應用要求,在用于實際量測數(shù)據(jù)分析時,同樣取得了很好的辨識效果,證明該方法能夠有效分解信號及準確提取低頻振蕩模態(tài)參數(shù)。
以類噪聲數(shù)據(jù)為基礎(chǔ),應用改進VMD算法進行數(shù)據(jù)的預處理,提出一種基于DFT曲線擬合的電力系統(tǒng)低頻振蕩信號識別方法。得出以下結(jié)論:
1)利用改進VMD 分解方法可有效消除直流分量或趨勢項,并能準確提取出低頻振蕩信號,抗噪性好;
2)利用自相關(guān)函數(shù)保持原信號振蕩模態(tài)參數(shù)特性,提出基于自相關(guān)函數(shù)DFT曲線擬合的模態(tài)參數(shù)辨識方法,通過DFT峰值個數(shù)可確定信號所含低頻振蕩模式的數(shù)量,運行時效性強;
3)采用模擬PMU 數(shù)據(jù)和某電網(wǎng)的實測數(shù)據(jù)計算分析,驗證了所提方法的準確性和有效性。該方法可用于環(huán)境激勵下的電力系統(tǒng)低頻振蕩在線模態(tài)參數(shù)辨識。