閆紅艷，Hwang Jin Kwon，高艷豐

（1. 河北工程大學水利水電學院，河北省 邯鄲市 056038；2. 韓國又石大學能源工程系，韓國鎮(zhèn) 川365-803）

0 引言

隨著電網遠距離、大容量輸電的實施，高放大倍數勵磁裝置的使用以及大規(guī)模新能源的接入，低頻振蕩時有發(fā)生，頻率范圍一般為0.1～2.5 Hz[1]，嚴重危害到系統正常運行。低頻振蕩模態(tài)識別是分析低頻振蕩產生原因的基礎，也是實現電力系統實時高效控制和風險預警的關鍵之一[2-5]。

基于相量測量單元(phasor measurement unit，PMU) 的電力系統廣域測量系統(wide area monitoring systems，WAMS)能夠實時地記錄電力系統的動態(tài)行為，按照量測信號的辨識方法可分為基于大擾動后自由振蕩響應信號的方法和基于環(huán)境激勵下隨機響應信號的方法[6-10]。大擾動數據包含強低頻振蕩分量，能夠準確識別區(qū)域間模式，但這類事件很少發(fā)生，很難進行連續(xù)的模態(tài)識別。相反，負荷隨機波動引起的小擾動激勵(又稱類噪聲數據)時刻存在、易于采集、數據豐富，可及時準確地反映當前系統的運行特性[11]。Prony 算法及其改進算法在分析大擾動激勵下的響應信號應用較廣泛，但Prony 算法在處理噪聲信號和具有時變特性非平穩(wěn)信號的能力不理想[12-13]。基于環(huán)境激勵下隨機響應信號的方法大多在信號的自回歸(auto regressive，AR)模型[10]或自回歸滑動平均(auto regressive moving average，ARMA)模型上發(fā)展起來的[14-16]。文獻[17-18]采用最小二乘法求取電力系統隨機響應信號AR 模型的最優(yōu)參數，進而獲得低頻振蕩的頻率和阻尼比信息，但是該算法存在數值不穩(wěn)定的問題。文獻[19]提出了用Yule-Walker 方程解AR 模型的功率譜分析方法；文獻[20]提出了一種基于數學形態(tài)學自回歸移動平均(mathematical morphology，MM-ARMA)算法的辨識方法。文獻[21]釆用Bootstrap 法確定YW 法辨識結果的置信區(qū)間。改進擴展的Yule-Walker(modified extended Yule-Walker，MEYW)方法[22]廣泛用于基于環(huán)境數據自相關函數的AR 模型中，但AR 模型的階數較難確定，存在2 個頻率相近的振蕩模式無法區(qū)分辨識的現象。

隨機子空間方法(stochastic subspace identification，SSI)從原理上既適用于自由振蕩信號分析，也適用于由環(huán)境激勵引起的系統隨機響應信號分析，抗噪性強，模態(tài)識別結果較準確，但計算速度較慢，很難反映信號的時變特性[7,23]。低頻振蕩產生的直流分量會降低模態(tài)識別的準確性，上述方法在使用過程中數據均需進行去除直流預處理。希爾伯特-黃變換(Hilbert-Huang transform，HHT)算法是分析非線性、非平穩(wěn)功率振蕩的常用工具，雖然能分解出直流分量，但存在端點效應現象[24]和本征模函數篩選效果不理想等固有缺點，辨識結果誤差率較大。文獻[25]提出集合經驗模態(tài)分解(ensemble empirical mode decomposition，EEMD)算法，但該算法存在計算復雜、耗時長等問題。

變分模態(tài)分解(variational mode decomposition，VMD)是一種新的非遞歸的信號分解方法，除了能濾除直流外，還能克服經驗模態(tài)分解(empirical mode decomposition，EMD)和EEMD 方法的缺陷，噪聲魯棒性較好。本文采用基于模態(tài)相關系數的VMD 進行類噪聲數據處理，直接去除中心頻率為0 的直流部分并分離出低頻振蕩信號；低頻振蕩信號模態(tài)參數辨識利用信號頻率差序列的自相關函數表示為與低頻振蕩函數形式一樣的指數衰減正弦函數的線性組合，該函數的固有頻率和阻尼比等于振蕩模態(tài)的固有頻率和阻尼比，用不同時間頻率數據之間的差序列減少測量噪聲的影響。通過在頻域中峰值頻率附近的離散傅里葉變換(discrete Fourier transform，DFT)進行曲線擬合來估計其拉普拉斯變換系數，然后計算得出模態(tài)分量的參數。通過與MEYW方法的識別結果相比較，驗證了所提識別方法的準確性。最后，采用基于ARMA 模型產生PMU 頻率數據和某實際系統的實測頻率數據驗證方法的可行性和有效性。

1 改進的變分模態(tài)分解方法

1.1 VMD原理

VMD 是Dragomiretskiy 提出的非遞歸信號分解方法，實質是變分問題，根據預設模態(tài)分量個數對信號進行分解[26]。該方法通過一個自適應維納濾波器組將原始信號f(x)分解為K個中心頻率為ωk的模態(tài)函數uk，其中K為預設模態(tài)分量個數。VMD算法在抗噪聲和非平穩(wěn)信號處理方面具有較好的性能和較高的運算效率，可以分解出直流分量。

為了得到具有一定帶寬頻率的K個模態(tài)分量，通常對每個模態(tài)函數uk進行Hilbert 變換得到邊際譜：

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

式中：{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；

1.2 基于相關系數的模態(tài)個數設定

VMD 本身不具有自適應性，模態(tài)分解數K值的設置是信號VMD 分解的關鍵環(huán)節(jié)，對分解效果影響較大。K設定值大于待分解信號所含的固有模態(tài)（intrinsic mode function，IMF）分量個數，則會在最終結果中引入虛假模態(tài)分量，影響對原始低頻振蕩信號的分析；相反，K設定值過小，則將導致信號分解不完全，即振蕩信號中含有的重要模態(tài)沒有被完全分解出來。

本文利用模態(tài)分量的相關系數來確定VMD分解個數的方法。首先取K值為2，計算各個分量之間的相關系數，判斷各分量之間是否存在頻率混疊現象，自適應確定模態(tài)個數。分量x1(n)和分量x2(n)的相關系數定義為

基于模態(tài)相關系數的VMD 算法具體步驟如下：

1）初始化K= 2，并用VMD 算法處理原始信號；

2）計算各個模態(tài)分量之間的相關系數，提取其中最大的相關系數；

3）K→K+1，并根據步驟2）更新模態(tài)之間的最大相關系數；

4）重復步驟3），直到最大相關系數超過閾值。經過大量實驗分析，最大相關系數的閾值選取0.1較為合適[27]。

根據文獻[26]，當懲罰參數α=2 000時，它可以滿足大多數工作條件的實際需求。試驗表明，改變α時影響很小，因此取默認值，當存在直流分量時，DC=1，這意味著可以分解出中心頻率為0 的分量。本文首要任務是要濾除直流分量，所以DC取1，其余參數取默認值。

2 電網區(qū)域間低頻振蕩模態(tài)辨識

2.1 類噪聲數據的ARMA模型

負荷的隨機變化即使很小也會影響到供需之間的不平衡，通常將這種干擾看作高斯白噪聲[28]，在電力系統類噪聲數據中用直流分量d(t)和低頻振蕩模式s(t)表示。同時，考慮到WAMS 數據中含有各種測量噪聲w(t)，因此，電力系統正常運行時，PMU數據的信號模型可表示為

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

式中v(t)=-u˙(t)，隨著v(t)的強度變強，由于電力系統中的調速器或頻率控制，隨著時間的推移，v(t)的強度可能會逐漸不同于其實際值。

擾動v(t)引起電力系統頻率中的低頻振蕩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]的功率比。在模擬產生數據和模態(tài)識別時，可將式(11)和式(13)的連續(xù)時間方程轉換為離散時間方程，用脈沖響應不變法[32]進行轉換，表示成ARMA模型為：

2.2 模態(tài)參數辨識

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

式中ck和ek分別為自回歸部分和滑動平均部分模型參數。

由信號的自相關特性可知，自相關分析能有效消除信號中的噪聲，且能保留原信號函數的頻率特征[33-34]，隨著時間的延長，噪聲信號自相關函數值將很快衰減至0，因此，可以采用自相關函數代替原函數進行模態(tài)辨識。

故sw[n]的自相關函數可以表示為

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

式中Re{·}代表復數的實數部分。

2.3 頻域的曲線擬合過程

公式(19)自相關函數rw[m]進行DFT的結果為

為了減少曲線擬合中的噪聲影響和縮短計算時間，希望在fk附近使用Rw[m]，fk是低頻振蕩模態(tài)功率最集中的地方。最好在以下頻率范圍上進行曲線擬合：|f^k-f|≤fc，k=1，2，…，K。

本文中將fc設為0.05 Hz，其寬度足以包含每個峰值的頻率范圍。每一個峰值曲線擬合采樣的最大整數小于2fc/△f。角頻率采樣點表示為ωk,m和ωk,(m+1)=ωk,m+2π△f。m=1,2,…,M是第K個模態(tài)峰值的曲線擬合。

信號VMD 分解后，取低頻振蕩頻率范圍內IMF分量，并對此分量做自相關計算，按式(24)中R(s)的系數ak和bk來辨識模式信號各參數的值。整個過程通過將R(s)擬合到頻域中的R[m]中來實現。R(s)在頻域中可以表示為

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

2.4 結果評價

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

式中：r[n]為測量信號；r^[n]為曲線擬合重構信號；SNRe的結果越大，表示擬合信號與原始信號擬合的效果越好。

本文采用VMD分解濾除原始數據中的直流分量，并提取出系統低頻振蕩信號，然后做自相關計算消除噪聲，低頻振蕩參數由自相關函數的拉普拉斯變換計算得到，而拉普拉斯變換式(22)的系數通過信號自相關函數DFT的曲線擬合來估計，最后由ARMA模型產生的類似真實PMU的頻率數據和實測PMU頻率數據驗證方法的有效性和準確性。

3 實例分析

3.1 模擬數據算例分析

由式(9)可知，PMU數據由3部分組成，為了驗證所提方法的可行性和有效性，利用式(15)和式(16)的ARMA模型模擬產生PMU頻率數據進行驗證。由ARMA模型可通過程序模擬產生離散的PMU數據，直流分量[27]取值如表1所示，2個低頻振蕩分量取值如表2 所示。模擬采樣時間間隔為0.02 s，加入不同信噪比的高斯白噪聲來模擬測量噪聲。

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

表2 低頻振蕩分量的參數Tab.2 Parameters of low-frequency oscillation components

模擬數據測量噪聲信噪比取30 dB 高斯白噪聲時，對模擬數據先經過低通濾波處理，低通濾波是基于高速采樣頻率50 Hz設計的2階巴特沃茲低通濾波器，截止頻率為5 Hz，濾除高頻分量可減少VMD 分解個數，增加運行速度。采用相關系數的VMD算法確定分量個數K，分量之間的最大相關系數如表3所示。由表3可知，當K取4、5時，信號經VMD算法分解之后分量之間的最大相關系數均大于閾值0.1；而當k取2、3時，最大相關系數皆小于閾值0.1。故取K值為3，VMD分解結果如圖1所示。

表3 不同K值的最大相關系數Tab.3 Maximum correlation coefficients of K different values

由圖1可知，信號經VMD分解后的各個IMF呈現比較規(guī)范，彼此間沒有模態(tài)混疊現象，各個頻段分離效果較好。其中IMF0是中心頻率為0的直流分量，其余的是不同頻率范圍的主導振蕩頻率，為后續(xù)準確辨識出低頻振蕩特征參數提取提供了理想的模態(tài)分量。

圖1 VMD模態(tài)分量及頻譜Fig.1 The spectrum and VMD modal component

與EEMD算法相比，VMD算法具有較好的優(yōu)越性，EEMD 分解結果如圖2 所示。從圖中可以看出，經過EEMD分解后得到13個IMF，右側為對應頻譜。分解得到的模態(tài)個數遠遠多于原始信號含有的振蕩分量個數，且耗時很長，同時出現了模態(tài)混疊，無法準確反映原始信號的低頻振蕩分量，嚴重影響參數提取的準確度。

以圖1 中IMF1 和IMF2 信號為例進行分析，對IMF1 分別利用頻差序列DFT 的曲線擬合法和MEYW法進行振蕩分量模態(tài)辨識，辨識算法采樣頻率為50 Hz，自相關函數的有效持續(xù)時間設置為20 s，數據窗長為5 min，相鄰數據窗間隔為1 min，數據總長為15 min。得到3～8 min 的頻率偏差波形和20 s頻率偏差的自相關擬合曲線如圖3所示，擬合DFT 幅值和角度如圖4 所示?？梢钥闯鰯M合曲線非常接近辨識信號的曲線，辨識結果準確度高。圖5給出IMF1分量的頻率和阻尼比的辨識結果，從圖中可以看出2 個方法的頻率偏差比較小，基于MEYW法的阻尼率值偏小，且受噪聲影響某個時間段誤差會比較大，而基于DFT的曲線擬合法變化比較平穩(wěn)，抗噪性好。將各滑動窗口內辨識得到的模態(tài)參數取平均值得到辨識結果。采用同樣的方法，將測量噪聲設成10 dB 高斯白噪聲時，K值為4，并對低頻振蕩頻率范圍內的分量進行模態(tài)辨識，模擬數據辨識結果如表4所示。

圖3 模擬數據IMF1分量的頻率偏差波形及自相關函數Fig.3 Frequency deviation waveform and autocorrelation function of IMF1

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

圖5 模擬數據IMF1頻率和阻尼率的辨識結果Fig.5 Frequencies and damping rates identification results of IMF1

表4 模擬數據的辨識結果Tab.4 Identification results of simulated data

由表4可以看出，考慮量測噪聲影響時，本文提出的辨識法能給出較為準確的頻率值，而阻尼比和振蕩幅值受系統運行方式的影響較大，所以信噪比越小，系統阻尼比和振蕩幅值波動程度相對越大，但本文提出的方法波動相對較小，說明所提方法準確度高和抗噪性好。因此，文中采用的基于DFT 曲線擬合的辨識法比目前應用廣泛的在線辨識MEYW法更準確，擬合精度SNRe如圖6所示。

對比圖6 中2 個分量的擬合精度SNRe可以看出，弱阻尼辨識結果比高阻尼辨識結果更準確，阻尼比越小，振蕩平息時間也越長，一段時間后該模式信號分量相對于噪聲仍占據主導地位，同時也可看出環(huán)境噪聲越弱，信噪比越高，辨識精度越高。

圖6 模擬數據擬合精度圖Fig.6 Fitting precision diagram

該方法的計算時間與VMD分解個數、采樣頻率、數序列長度、數據窗滑動步長以及時間窗長度直接相關，VMD計算時間占比較大，但相對于EEMD 計算時間能減少很多，同時低通濾波后可減少分解個數，提高運行速度。由于MEYW法中AR模型的階數比DFT擬合法的階數高，相比較，本文辨識方法運行時間短。以MATLAB 2018 版進行編程，完成上述PMU模擬數據分析，SNR為30 dB 時低通濾波后VMD 分解個數為3，計算時間為6.372 7 s，基于DFT 曲線擬合辨識運行時間為0.927 8 s，耗時7.518 9 s，占空比為0.835 4%，完全滿足在線應用要求。

3.2 實測數據模態(tài)辨識

以某實際系統PMU錄波數據為例，分析系統低頻振蕩特征。實際系統采樣頻率為0.033 3 s，選取時長為15 min 的部分量測信號，實測數據減去額定頻率60 Hz，得到實際頻率的波動如圖7所示。

圖7 某電網實測頻率波動Fig.7 The measured frequency fluctuation of a power grid

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

表5 實測數據辨識結果Tab.5 Identification results of measured data

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

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

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

圖11 實測數據擬合精度Fig.11 Fitting accuracy of measured data

在實測數據辨識過程中，采樣頻率較低，為30 Hz，同樣在MATLAB 2018 版本上運行程序。低通濾波后VMD 分解個數為4，運行時間為8.404 8 s，總耗時為9.551 s，占空比為1.061%，提高了機電小干擾穩(wěn)定評估的實時性。

以上結果表明，由于VMD 算法和基于DFT的曲線擬合法都具有較好的噪聲魯棒性，參數辨識結果精度較高，辨識曲線平滑，符合實際的波動情況，并且計算速度滿足在線應用要求，在用于實際量測數據分析時，同樣取得了很好的辨識效果，證明該方法能夠有效分解信號及準確提取低頻振蕩模態(tài)參數。

4 結論

以類噪聲數據為基礎，應用改進VMD算法進行數據的預處理，提出一種基于DFT曲線擬合的電力系統低頻振蕩信號識別方法。得出以下結論：

1）利用改進VMD 分解方法可有效消除直流分量或趨勢項，并能準確提取出低頻振蕩信號，抗噪性好；

2）利用自相關函數保持原信號振蕩模態(tài)參數特性，提出基于自相關函數DFT曲線擬合的模態(tài)參數辨識方法，通過DFT峰值個數可確定信號所含低頻振蕩模式的數量，運行時效性強；

3）采用模擬PMU 數據和某電網的實測數據計算分析，驗證了所提方法的準確性和有效性。該方法可用于環(huán)境激勵下的電力系統低頻振蕩在線模態(tài)參數辨識。