項 頌, 蘇 鵬, 吳曉丹, 王 楊, 蔣小龍

(1. 國網(wǎng)內(nèi)蒙古東部電力有限公司, 內(nèi)蒙古 呼和浩特 010010; 2. 四川大學(xué)電氣工程學(xué)院, 四川 成都 610065)

1 引言

近年來,以風(fēng)電和光伏為代表的大量可再生能源通過電力電子設(shè)備接入電網(wǎng),逐步形成了以高比例可再生能源和高比例電力電子設(shè)備為特征的“雙高”電力系統(tǒng)[1,2]。而電力電子設(shè)備帶來的弱抗擾、低慣性以及強耦合等特點已嚴重影響了電力系統(tǒng)的安全穩(wěn)定運行,導(dǎo)致寬頻振蕩事故頻發(fā),如美國德州電網(wǎng)的雙饋風(fēng)電場發(fā)生了多起20 Hz次同步振蕩(SubSynchronous Oscillation, SSO)事故[3];我國河北沽源地區(qū)出現(xiàn)了多次3～10 Hz次同步振蕩現(xiàn)象[4];法國-西班牙聯(lián)網(wǎng)工程發(fā)生了1.4 kHz的振蕩事故;我國新疆哈密地區(qū)直驅(qū)風(fēng)電機組發(fā)生了20～80 Hz的次/超同步振蕩事故[5];云南魯西柔直工程發(fā)生了1.27 kHz的高頻振蕩事故。若振蕩事故沒有得到有效抑制,可能造成設(shè)備損毀、可再生能源機組脫網(wǎng),甚至傳遞至用戶側(cè)造成大范圍的停電事故。隨著中國“雙碳”目標的提出,可再生能源滲透率將進一步提高,電網(wǎng)抵御寬頻振蕩的能力將進一步受到挑戰(zhàn)。

目前,國內(nèi)外針對電力系統(tǒng)寬頻振蕩的分析方法以建模分析為主,如阻抗建模法和狀態(tài)空間法。盡管該類方法具有堅實的理論基礎(chǔ),但其分析精度在實際應(yīng)用中極大地受制于系統(tǒng)模型的準確性?；诒O(jiān)測數(shù)據(jù),分析人員可以從全局的角度進行振蕩風(fēng)險評估、振蕩事故反演、振蕩源定位等[6-8]。文獻[9]提出,在系統(tǒng)模型不明確的情況下,構(gòu)建寬頻振蕩的廣域監(jiān)測系統(tǒng)是分析寬頻振蕩產(chǎn)生機理與傳播特性的有效手段。同時,構(gòu)建寬頻振蕩監(jiān)測系統(tǒng)也是實現(xiàn)振蕩實時預(yù)警、控制、保護的硬件基礎(chǔ)。

當前的廣域測量系統(tǒng)以相量測量單元(Phasor Measurement Unit, PMU)為基礎(chǔ),只關(guān)注工頻分量[9],難以滿足寬頻振蕩監(jiān)測需求,為此,可以利用子站PMU的波形記錄功能或改進廣域監(jiān)測系統(tǒng)(Wide-Area Measurement System, WAMS)動態(tài)監(jiān)測算法將同步測量能力擴大到寬頻范圍。目前寬頻振蕩監(jiān)測的常用方法包括傅里葉變換算法、小波變換算法、希爾伯特黃變換算法(Hilbert-Huang Transform, HHT)、卡爾曼濾波算法、Prony算法等。文獻[10-12]均采用加窗的插值傅里葉變換(Interpolated Discrete Fourier Transform, IpDFT)結(jié)合濾波器進行寬頻振蕩參數(shù)估計,對穩(wěn)態(tài)信號的估計比較準確,但是難以處理非平穩(wěn)信號。文獻[13,14]采用了小波算法處理時變信號,但是在遇到多頻率分量時存在分辨不高的問題。文獻[15,16]使用了HHT算法,但是并沒有完全解決端點效應(yīng)及模態(tài)混疊的問題。文獻[17-19]利用卡爾曼濾波算法對諧波進行跟蹤,但是所需參數(shù)較多,計算量較大。文獻[20,21]將形態(tài)濾波和小波算法及Prony算法結(jié)合,提高了一定的抗噪性能,但沒有解決模型階數(shù)難以確定及計算量較大的問題。除此之外,上述算法在處理強時變特性的信號時均會出現(xiàn)較大誤差,甚至?xí)耆?因而無法直接應(yīng)用于復(fù)雜多變的振蕩信號中。

泰勒傅里葉變換(Taylor-Fourier Transform, TFT)算法因采用動態(tài)相量模型而具有更為優(yōu)越的動態(tài)性能[22],因此本文將TFT算法用于寬頻振蕩監(jiān)測中,但TFT算法易受到噪聲干擾,常需要增加數(shù)據(jù)窗長[23]或者進行濾波預(yù)處理[24]。經(jīng)作者大量測試得知,常規(guī)濾波算法難以解決TFT算法的抗噪問題或以大幅提升計算時間為代價,部分濾波算法甚至?xí)档蚑FT算法的精度,因此增加窗長是更為有效的抗噪方案。而TFT算法在噪聲較小時并不需要較長的窗長,故本文提出根據(jù)噪聲強度自適應(yīng)調(diào)整窗長的方案。徑向基(Radial Basis Function, RBF)神經(jīng)網(wǎng)絡(luò)具有較強輸入和輸出的映射能力,可以高精度逼近任意非線性函數(shù),且有拓撲結(jié)構(gòu)緊湊、收斂速度快等優(yōu)點,因此本文將其應(yīng)用于噪聲強度估計。最后,本文通過大量仿真測試及河北沽源和新疆哈密的實測數(shù)據(jù)進行了驗證,結(jié)果表明文中方法具有較好的可行性和準確性。

2 泰勒傅里葉變換

電力系統(tǒng)寬頻振蕩信號的時域模型為:

(1)

式中,Am(t)、fm、φm(t)分別為第m個信號分量的峰值、頻率和初相位;N(t)為系統(tǒng)的噪聲;M為信號分量總數(shù);t為時間。

以頻率fs對x(t)進行采樣,經(jīng)采樣后的離散相量模型可表示為:

(2)

pm(n)=Am(nTs)ejφm(nTs)

(3)

式中,ωm為信號的角頻率;pm(n)為第m個信號分量的相量,n=-Nw,…,0,…,Nw;Ts為采樣周期。

將相量pm(n)在零點進行k階泰勒級數(shù)展開:

(4)

式中,k為泰勒展開的階數(shù)。

根據(jù)歐拉公式可以將式(2)寫成矩陣形式:

(5)

式中,X為寬頻振蕩信號的采樣序列;B為泰勒傅里葉的基向量矩陣;P為相量的各階導(dǎo)數(shù)及其共軛項所組成的列向量,表達如下:

X=[x(-Nw) …x(0) …x(Nw)]T

(6)

P=[ph1…phm…phM]T

(7)

B=[B1…Bm…BM]

(8)

phm=(pm(k)(0),…,pm(0),pm(0)*,…,pm(k)(0)*)

(9)

(10)

式中,*為共軛運算符。

根據(jù)矩陣理論,利用最小二乘法便可求解P。

(11)

式中,H為Hermitian運算符;W為權(quán)重對角矩陣,用于降低高頻干擾的影響。寬頻振蕩信號的峰值、相位和頻率計算為:

(12)

式中,fm0為第m個分量的頻率初始值;abs(·)為取模運算符;angle(·)為求角運算符。

3 TFT算法的影響因素

3.1 初始頻率

由基向量矩陣B的構(gòu)成可知,TFT算法極其依賴初始代入頻率,錯判頻率個數(shù)甚至?xí)?dǎo)致TFT算法直接失效。相反,如果代入頻率接近實際頻率,TFT算法的辨識精度就會得到保證。因此,本文采用了離散傅里葉變換(Discrete Fourier Transform, DFT)算法來確定初始頻率,為了降低頻譜泄露以及柵欄效應(yīng)的影響,使用了窗函數(shù)及插值算法,窗函數(shù)選擇了通帶平穩(wěn)、阻帶衰減快的Hanning窗,插值算法選擇了三譜線插值算法[25]。

考慮到加窗插值算法無法完全消除柵欄效應(yīng)和頻譜泄露的影響,初始頻率的準確性仍可能較低,因此將由式(9)得到的頻率回代到式(10)中迭代更新B矩陣并重新用式(9)計算相應(yīng)參數(shù)。理論上來說,迭代的次數(shù)越多,TFT算法越精確,但是計算效率也會相應(yīng)降低。

3.2 窗口長度

TFT算法在計算向量P時所采用的加權(quán)最小二乘法在窗口較長時計算量大、運算時間長,在實際應(yīng)用時需根據(jù)具體情況盡可能地縮短所需數(shù)據(jù)窗長。然而向量P中包含了較多的高階導(dǎo)數(shù)項,高階導(dǎo)數(shù)項對噪聲比較敏感,在噪聲較大時需要足夠多的數(shù)據(jù)才能保證其準確性,因此需要同時考慮噪聲強度及計算量來確定TFT算法的數(shù)據(jù)窗長,故提出了基于噪聲強度的窗長自適應(yīng)方法,步驟如下:

(1)利用RBF神經(jīng)網(wǎng)絡(luò)算法擬合信號。

(2)將采樣波形與擬合波形做序列減法,得到噪聲序列,然后估計信噪比。

(3)根據(jù)信噪比計算數(shù)據(jù)窗長。

3.2.1 擬合去噪

人工神經(jīng)網(wǎng)絡(luò)相比傳統(tǒng)擬合算法具有對系統(tǒng)模型的低依賴性、強大自適應(yīng)能力以及容錯性高等優(yōu)點,且非常適合處理復(fù)雜的非線性問題,可以逼近任意函數(shù),因而本文采用神經(jīng)網(wǎng)絡(luò)算法對信號進行擬合。在各類神經(jīng)網(wǎng)絡(luò)算法中,RBF神經(jīng)網(wǎng)絡(luò)具有結(jié)構(gòu)簡單、速度快且函數(shù)擬合精度高的優(yōu)點,非常適用于擬合去噪[26]。

RBF神經(jīng)網(wǎng)絡(luò)是一種三層前向網(wǎng)絡(luò),由輸入層、隱含層及輸出層組成,結(jié)構(gòu)如圖1所示。

圖1 RBF神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.1 Structure of RBF neural network

圖1中,X為輸入層;xi為輸入向量;R為隱含層;ri為神經(jīng)元;wj為權(quán)重向量;F為輸出層空間。RBF神經(jīng)網(wǎng)絡(luò)采用徑向基函數(shù)作為激活函數(shù),從而直接將輸入向量空間映射到隱含層,因此隱含層是非線性的,輸出層則為線性的。

(13)

式中,f(x)為輸出層的輸出,由隱含層神經(jīng)元輸出的加權(quán)求和可得;φ(‖x-xj‖)為隱含層的激活函數(shù),多為高斯函數(shù),如下所示:

(14)

式中,σ為方差;x0為徑向基函數(shù)的中心點。

令輸出層的輸出為輸入向量的期望響應(yīng),記為:

FT=[f1f2…fu]

(15)

同時,令:

φij=φ(‖xi-xj‖)

WT=[w1w2…wv]

(16)

則式(13)可以改寫為:

(17)

即:

ΦW=F

(18)

如果Φ可逆,則可得:

W=Φ-1F

(19)

本文采用正則化RBF網(wǎng)絡(luò),其隱含層節(jié)點數(shù)與輸入樣本數(shù)相同,且各輸入樣本均為徑向基函數(shù)的中心點,因此無需再另外設(shè)計算法。方差σ計算為:

(20)

式中,dmax為中心點之間的最大距離;b為隱含層節(jié)點數(shù)。

由于RBF神經(jīng)網(wǎng)絡(luò)會帶來運算量的增加,降低實時性,因此用于擬合的窗不宜過長,文中采用1個基波周期的數(shù)據(jù)用于噪聲估計。

3.2.2 噪聲強度估計

將原始波形與經(jīng)RBF神經(jīng)網(wǎng)絡(luò)擬合的波形做差即可得到噪聲數(shù)據(jù),信噪比便可根據(jù)式(21)得到,單位為分貝(dB)。

(21)

式中,PS和PN分別為信號和噪聲的功率。

3.2.3 確定窗長

信噪比越高,所需數(shù)據(jù)越少。對于采樣率為12 800 Hz的系統(tǒng),通過測試發(fā)現(xiàn),當信噪比為30 dB及以上時,TFT算法只需要5個基波周期即可保證高精度的參數(shù)估計。當信噪比較低時,窗長應(yīng)根據(jù)實際情況適當增加,但仍不建議超過10個基波周期以保證計算效率和動態(tài)性能。因此,窗長Win采用分段處理(遇小數(shù)則四舍五入):

(22)

文中所需的窗長一直保持在5～10個基波周期,在信噪比較低時,該方法會在一定程度上降低算法的動態(tài)實時性,但是削弱了噪聲對參數(shù)辨識結(jié)果的影響。因此,本文方法兼顧了動態(tài)實時性和準確性,在實際應(yīng)用中,窗長函數(shù)可根據(jù)應(yīng)用場景和需求進行適當調(diào)整。

3.3 算法實現(xiàn)流程

本文算法的整體流程如圖2所示,首先確定泰勒展開的階數(shù)以及TFT的迭代次數(shù),然后對采樣序列進行快速估計和噪聲強度估計分別得到初始頻率和窗口長度,最后基于TFT算法通過多次迭代計算出頻率、幅值和相位等參數(shù)。

圖2 算法流程圖Fig.2 Flowchart of algorithm

4 仿真結(jié)果及分析

4.1 性能評價標準

本文采用IEEE Std. C37.118.1-2011[27]標準中以評估相量測量算法性能的總體相量誤差(Total Vector Error, TVE)和頻率誤差(Frequency Error, FE)作為評價標準,TVE和FE的定義如下:

(23)

(24)

4.2 仿真信號測試

以次同步振蕩為例(基波存在頻率偏移)驗證本文算法的效果,仿真信號如式(25)所示,并加入了30 dB的高斯白噪聲,采樣率為12 800 Hz,泰勒傅里葉階數(shù)為2階,共迭代3次,波形如圖3所示。

圖3 含次同步振蕩波形圖Fig.3 Waveform plot with subsynchronous oscillation

(25)

30 dB的信噪比對應(yīng)的窗長為5個基波周期,擬合結(jié)果如圖4所示,工頻分量的總體相量誤差百分比和頻率誤差分別為0.16%和0.003 Hz,SSO的則為1.77%和0.33 Hz。

圖4 擬合對比圖Fig.4 Comparison between waveforms

為進一步驗證本文方法在復(fù)雜工況下的效果,本文參考了IEEE Std. C37.118.1—2011標準中的PMU測試標準,設(shè)計了如下仿真算例,且各仿真的噪聲水平均從15～60 dB變化。

(1) 幅值動態(tài)

在真實情況下,SSO發(fā)生的過程往往比較迅猛,本例以衰減指數(shù)形式模擬振蕩時的幅值動態(tài)變化,衰減阻尼為1,信號模型為式(26),仿真結(jié)果如圖5所示。

(26)

(2) 頻率斜坡動態(tài)

在實際的振蕩事故中,SSO頻率可能具有時變性,故以斜坡函數(shù)為例進行仿真,信號模型為式(27),SSO分量的頻率分量疊加了線性項10t,仿真結(jié)果如圖6所示。

圖6 頻率斜坡動態(tài)下的誤差Fig.6 Errors under frequency ramp test

(27)

(3) 頻率調(diào)制動態(tài)

在實際電力系統(tǒng)中,SSO頻率的變化并非一定為線性變化,此處對頻率進行正弦調(diào)制,調(diào)制信號見式(28),頻率疊加了正弦項-10sin(2π×0.1t),仿真結(jié)果如圖7所示。

(28)

(4) 多模式寬頻振蕩

寬頻振蕩中可能存在多個信號分量,比如會有諧波和間諧波同時存在的情況,而且各信號分量的動態(tài)過程都各有不同,此時的工況愈發(fā)復(fù)雜,對監(jiān)測算法的要求也更高。本例的仿真信號為式(29),具體參數(shù)見表1,仿真結(jié)果如圖8所示。

表1 信號參數(shù)Tab.1 Parameters of signals

圖8 多模式寬頻振蕩下的誤差Fig.8 Error under multi-mode broadband oscillation

(29)

4.3 算法性能討論

(1)從仿真結(jié)果中可以看出,本文算法的動態(tài)性能優(yōu)越,幾乎可以應(yīng)對任何動態(tài)變化過程,且在噪聲較強時仍具有較高準確性。因TFT算法對泰勒展開項進行了二階截斷,所以當振蕩動態(tài)變化較快時誤差也會有所增大。在實際應(yīng)用時,如果算力允許,可增加泰勒展開項階數(shù)以進一步增強TFT算法的動態(tài)性能。

(2)工頻分量的相量誤差和頻率誤差遠小于其他信號分量,這是因為工頻分量的幅值較大,受噪聲影響較小。根據(jù)文獻[27]中對工頻相量的測量標準(最大總體相量誤差為1%,最大的頻率誤差為0.01 Hz),由圖5～圖8可知,本文算法滿足要求。

(3)TFT算法極其依賴初始頻率,而振蕩信號在振蕩初期可能因噪聲較大而無法從插值傅里葉算法中識別出來,這有可能導(dǎo)致TFT算法崩潰,因此在使用本算法時可對采樣信號進行預(yù)處理。

(4)數(shù)據(jù)窗長對TFT算法也具有較大的影響,在信噪比較低時,若采樣率較低,可不限于5個基波周期和10個基波周期的最低限制和最高限制。若需提升動態(tài)實時性,也可犧牲準確性降低窗長自適應(yīng)范圍,具體可根據(jù)需要做取舍。

(5)需要指出的是,當測量低頻振蕩時,本文方法與其他基于傅里葉變換的方法相似[10],窗長須有所增加才能保證其準確性。本文當前采用的窗長方案僅適合分析5 Hz及以上的信號。經(jīng)試驗,若需分析更為低頻的信號,窗長需要適當增加,如15個基波周期才能達到較高的精確性。而高頻振蕩并不會影響本文算法的準確度。由于篇幅限制,分析結(jié)果不再展示。

(6)本文的測試環(huán)境為Matlab R2020b,硬件配置AMD Ryzen 5-4500U CPU@3.00 GHz,運行內(nèi)存為8 GB DDR4。當信號包含少量諧波或間諧波信號時,本文算法的計算效率較高,比如信號中僅包含次同步振蕩分量時,算法運行時間約為0.02 s。但是當信號存在大量的諧波或間諧波時,如本文多模式寬頻振蕩,此時運行時間約為0.8 s。其主要原因在于信號較多時,TFT算法中的矩陣B的維度會大幅增加,所以導(dǎo)致求解向量P會十分耗時,且RBF算法的擬合時間約為0.01 s?？紤]到本文所用的測試環(huán)境計算能力較弱,在實際應(yīng)用時,可選擇高性能硬件處理器進一步減少計算時間[24]。

4.4 算法對比

為兼顧辨識算法的穩(wěn)定性和準確性,目前工程實際應(yīng)用中多采用基于傅里葉變換的算法實現(xiàn)寬頻量測[10]。文獻[10]采用了漢寧窗和雙譜線插值算法解決傳統(tǒng)傅里葉變換中易出現(xiàn)的頻譜泄露和柵欄效應(yīng)問題。以頻率調(diào)制動態(tài)為例,傅里葉算法窗長固定為0.2 s,仿真設(shè)置與式(28)一致,對比結(jié)果如圖9所示。

圖9 算法對比Fig.9 Comparison with other algorithms

從圖9對比結(jié)果中可以看出,在辨識頻率非線性變化的次同步分量時,傅里葉算法辨識結(jié)果的誤差較大,這是因為傅里葉算法會對窗口內(nèi)的特征參數(shù)進行“平均”,所以該算法更適用于線性變化的信號。而對非線性信號的適應(yīng)性極差,在辨識穩(wěn)態(tài)工頻分量時,其辨識結(jié)果也不如本文算法,主要是因為傅里葉變換受噪聲和數(shù)據(jù)窗長的影響更大,而且還存在頻譜泄露及柵欄效應(yīng)等問題。

為展示RBF神經(jīng)網(wǎng)絡(luò)的擬合效果,本文將局部加權(quán)回歸平滑算法(Locally Weighted Regression Smoothing, LWRS)、切比雪夫Ⅱ型(ChebyshevII)數(shù)字濾波器以及小波算法與RBF神經(jīng)網(wǎng)絡(luò)進行比較,信號模型為:

(30)

其中,仿真信號添加了15 dB的高斯白噪聲。另外為方便觀察,采樣率降低為1 600 Hz。LWRS算法的擬合階數(shù)為2階、窗口長度為15;切比雪夫Ⅱ型濾波器的通帶頻率為80 Hz、通帶波紋為0.1 dB、阻帶頻率為81 Hz、阻帶衰減為80 dB;小波算法采用了dmey小波,分解層數(shù)設(shè)為3。擬合誤差對比如圖10所示。

從擬合結(jié)果和誤差對比圖可以看出,RBF神經(jīng)網(wǎng)絡(luò)的擬合誤差更小,信噪比估算更準確,而且沒有復(fù)雜的參數(shù)設(shè)置過程。綜上,在實際應(yīng)用中,本文算法具有一定的優(yōu)越性和實用性。

5 實測數(shù)據(jù)驗證

本節(jié)分別以河北沽源和新疆哈密地區(qū)真實發(fā)生的次同步及次/超同步振蕩案例對本文方法進行驗證。

沽源和哈密系統(tǒng)實測的電流分別如圖11和圖12所示,辨識結(jié)果分別如圖13和圖14所示。根據(jù)RBF擬合結(jié)果可知,沽源及哈密的噪聲強度分別約為30 dB和55 dB,但由于沽源與哈密的故障錄波儀的采樣率較低,分別只有3 200 Hz和1 200 Hz,故相應(yīng)增大了數(shù)據(jù)窗長(10周期)。

圖12 哈密某場站電流波形Fig.12 Current waveform of a wind farm in Hami

從辨識結(jié)果中能夠得知,沽源只發(fā)生了次同步振蕩事故,而哈密則同時存在次/超同步分量,而且從頻率和幅值變化的過程可以判斷出這是一對同源的次超同步分量。根據(jù)事故后的分析可知,沽源和哈密振蕩事故的發(fā)生原理不同,前者為雙饋風(fēng)機與固定串補之間的交互作用,而后者為永磁直驅(qū)風(fēng)機與弱電網(wǎng)之間的交互作用,但無論是哪種振蕩事故,短時間內(nèi)振蕩的幅值都會迅速增大,進而可能危害到系統(tǒng)的安全穩(wěn)定運行。

6 結(jié)論

本文提出了一種基于波形數(shù)據(jù)的寬頻振蕩監(jiān)測方法。首先,利用DFT算法得到初步的頻率參數(shù),并利用RBF神經(jīng)網(wǎng)絡(luò)進行擬合去噪,從而得到噪聲強度并自適應(yīng)確定數(shù)據(jù)窗長,根據(jù)所得窗長和初始頻率利用TFT算法迭代提取出相關(guān)振蕩參數(shù)。文中通過動態(tài)仿真測試和實際的振蕩事故案例進行了驗證,從測試結(jié)果可以看出本文采用的RBF神經(jīng)網(wǎng)絡(luò)算法具有較好的擬合效果,且窗長的自適應(yīng)處理可以保證TFT算法的計算效率和實時性。從參數(shù)辨識結(jié)果看本文方法的動態(tài)性能優(yōu)越,在各種工況下均有較強的適應(yīng)性。

本文方法不局限于振蕩監(jiān)測,還可根據(jù)實際需求拓展至其他領(lǐng)域,如諧波檢測、電能質(zhì)量等,并在未來實現(xiàn)基于廣域監(jiān)測系統(tǒng)的寬頻監(jiān)測、穩(wěn)定分析、風(fēng)險預(yù)警等。RBF神經(jīng)網(wǎng)絡(luò)已廣泛應(yīng)用于各工程領(lǐng)域中[28,29],本文所采用的基于RBF神經(jīng)網(wǎng)絡(luò)的擬合方法所帶來的計算負擔完全在可接受范圍內(nèi),但在振蕩分量較多且噪聲較大時,TFT算法的計算效率會受到影響,因此在實時監(jiān)測方面的應(yīng)用仍需進一步研究。