林偉斌，季天瑤，張祿亮

(華南理工大學(xué) 電力學(xué)院，廣東 廣州 510641)

目前，電網(wǎng)規(guī)模不斷擴(kuò)大，系統(tǒng)間的互聯(lián)程度不斷增大，越來越多的高倍率快速勵磁裝置投入使用，導(dǎo)致系統(tǒng)弱阻尼低頻振蕩出現(xiàn)概率更高[1-2]，極大地危害電網(wǎng)的穩(wěn)定性，也限制了互聯(lián)系統(tǒng)傳輸容量最大化。因此，監(jiān)控和分析低頻振蕩模態(tài)參數(shù)，對確保電力系統(tǒng)的安全穩(wěn)定具有重要意義。

低頻振蕩模態(tài)識別方法可以分為基于模型和基于量測信號方法[3]。由于電網(wǎng)結(jié)構(gòu)越來越復(fù)雜，基于模型的方法計算復(fù)雜，建模困難，會導(dǎo)致出現(xiàn)維數(shù)災(zāi)等問題，且無法實時計算參數(shù)，在實際應(yīng)用中受到限制。近年來，廣域測量系統(tǒng)已在電力系統(tǒng)中大規(guī)模使用，為基于量測信號的方法提供了數(shù)據(jù)支持[4]?；诹繙y信號的方法無需知道系統(tǒng)結(jié)構(gòu)，可以在線快速準(zhǔn)確地識別低頻振蕩參數(shù)，是一種更簡單直接的方法。隨著電網(wǎng)結(jié)構(gòu)愈加復(fù)雜，基于量測信號的方法逐漸成為低頻振蕩模態(tài)識別研究的重點。

根據(jù)激勵源的不同，量測信號可分為：由明顯擾動(如短路故障、投切大負(fù)荷)引起的暫態(tài)振蕩信號和正常運行時負(fù)荷隨機(jī)波動引起的類噪聲信號[5]?；跁簯B(tài)振蕩信號的方法可以準(zhǔn)確地提取所需的低頻振蕩模態(tài)參數(shù)，許多專家也進(jìn)行了相關(guān)研究并取得了良好的結(jié)果。采用暫態(tài)振蕩信號進(jìn)行低頻振蕩模態(tài)估計最具有代表性的方法是Prony算法，其算法簡便，識別結(jié)果令人滿意，但抗噪性能差，在含噪情況下識別精度差。文獻(xiàn)[6]使用了時域和頻域混合的方法進(jìn)行模態(tài)識別，該方法與Prony算法相比有更好的抗噪性能和更高的識別精度。此外，基于暫態(tài)振蕩信號的常用方法包括傅里葉變換[7]、小波變換[8]、希爾伯特-黃變換[9]。雖然基于暫態(tài)振蕩信號進(jìn)行低頻振蕩模態(tài)辨識效果良好，但暫態(tài)振蕩信號發(fā)生概率低，難以實時獲得，無法實時反映系統(tǒng)狀態(tài)，難以判斷正常運行期間的小干擾穩(wěn)定性。相反，由負(fù)荷波動引起的類噪聲信號幾乎一直存在，可實時獲取?；陬愒肼曅盘柕姆椒梢詫崟r監(jiān)控系統(tǒng)狀態(tài)，在發(fā)生明顯振蕩前識別出低頻振蕩的模態(tài)參數(shù)，從而實現(xiàn)對弱阻尼或負(fù)阻尼的預(yù)警，這是電力系統(tǒng)調(diào)度運行人員更為關(guān)心的，因此該方法在低頻振蕩的監(jiān)測中有更好的應(yīng)用前景。

近年來，基于類噪聲信號的低頻振蕩方法大量涌現(xiàn)，基于隨機(jī)子空間識別(stochastic subspace identification，SSI)的方法是其中的代表。該方法應(yīng)用系統(tǒng)的狀態(tài)方程進(jìn)行參數(shù)估計，能方便直接地從數(shù)據(jù)中獲取狀態(tài)，算法參數(shù)中僅系統(tǒng)的階次較難確定。為解決模型定階問題，文獻(xiàn)[10]采用穩(wěn)定圖法自動定階，文獻(xiàn)[11]采用經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition，EMD)與SSI相結(jié)合的方法進(jìn)行模態(tài)辨識，文獻(xiàn)[12]引入一致性指示值判別真實模態(tài)和虛假模態(tài)。但上述研究均是基于單一維度下的漢克爾矩陣得到的模態(tài)穩(wěn)定圖進(jìn)行模型定階。本文提出一種新型的基于SSI的低頻振蕩模態(tài)辨識方法。該方法是在SSI算法的基礎(chǔ)上，構(gòu)造2個不同維度的漢克爾矩陣進(jìn)行模態(tài)辨識，基于2個不同維度的漢克爾矩陣下得到的同階極點進(jìn)行匹配而得到超清穩(wěn)定圖，最后再對超清穩(wěn)定圖進(jìn)行系統(tǒng)聚類，獲得最終的真實模態(tài)參數(shù)。

本文方法不需人工進(jìn)行定階，相比于基于單一維度下的SSI方法能更有效地區(qū)分真實模態(tài)和虛假模態(tài)。采用傳遞函數(shù)構(gòu)造的合成數(shù)據(jù)以及3機(jī)9節(jié)點系統(tǒng)仿真數(shù)據(jù)對本文方法的有效性進(jìn)行驗證，并與傳統(tǒng)SSI算法進(jìn)行對比。

1 基于協(xié)方差的SSI算法

由環(huán)境激勵的自由度離散系統(tǒng)的狀態(tài)方程為[13]：

(1)

式中：A∈Rn×n為系統(tǒng)的狀態(tài)矩陣；C∈Rl×n為系統(tǒng)的輸入矩陣；xk∈Rn和yk∈Rl分別為系統(tǒng)k時刻的狀態(tài)量和輸出量；wk∈Rn是驅(qū)動過程中的白噪聲擾動量；vk∈Rl為測量噪聲；wk和vk均為不可測量的信號，假設(shè)它們?yōu)榛ゲ幌嚓P(guān)、均值為零的白噪聲平穩(wěn)序列；n為系統(tǒng)階數(shù)，在本文中為低頻振蕩模態(tài)數(shù)；l為系統(tǒng)輸出數(shù)量。

在離散域中，對狀態(tài)矩陣A進(jìn)行特征值分解：

A=ψΛψ-1.

(2)

式中：Λ=diag(λs)，s=1,2,…,n，λs為分解得到的第s個離散模態(tài)特征值；ψ為特征矩陣。對應(yīng)的第s個連續(xù)時間特征值λs,c及低頻振蕩模態(tài)參數(shù)為：

(3)

(4)

式中：Δt為采樣時間間隔；fs為低頻振蕩模態(tài)頻率；ξs為低頻振蕩阻尼比；ωs為低頻振蕩模態(tài)角速度；Re表示取特征值實部。

SSI是一種用于模態(tài)參數(shù)估計的高效系統(tǒng)識別工具，具有數(shù)值簡單性和魯棒性的優(yōu)點。根據(jù)投影矩陣的不同，SSI分為DATA-SSI和協(xié)方差驅(qū)動的SSI[14]。在保留原有信號的信息的情況下，DATA-SSI采用QR分解對數(shù)據(jù)量進(jìn)行壓縮，而協(xié)方差驅(qū)動的SSI采用托普利茨矩陣對數(shù)據(jù)量進(jìn)行壓縮。在實際運算中，協(xié)方差驅(qū)動的SSI可采用快速傅里葉變換進(jìn)行數(shù)據(jù)壓縮，而DATA-SSI采用的是QR分解，運算速度較慢。為了滿足快速計算和實時估計低頻振蕩參數(shù)的要求，本文采用協(xié)方差驅(qū)動的SSI來進(jìn)行模態(tài)參數(shù)估計，其數(shù)學(xué)過程如下[15]。

首先由系統(tǒng)輸出量構(gòu)造漢克爾矩陣：

(5)

(6)

(7)

式中a、b分別為漢克爾矩陣的行數(shù)和列數(shù)。理論上，為滿足統(tǒng)計分析的需求，b→，即獲得的系統(tǒng)輸出量數(shù)據(jù)越多時，識別結(jié)果越準(zhǔn)確。但實際中無法做到數(shù)據(jù)量無窮大，且數(shù)據(jù)越多計算負(fù)擔(dān)越大。相關(guān)文獻(xiàn)表明，當(dāng)a> round (n/l)(round(·)代表向上取整)，且b>20a時，SSI算法的識別結(jié)果能滿足精度要求。在本文中，為了在計算精度和計算成本之間取得一個較好的平衡，b的取值對應(yīng)時間長度為10 min的數(shù)據(jù)。下標(biāo)“p”表示“過去”，“f”和“f，2”表示“未來”，即Yp為“過去”漢克爾矩陣，Yf和Yf,2為第1和第2個“未來”矩陣。

根據(jù)式(5)—(7)構(gòu)造托普利茨矩陣：

(8)

(9)

對式(8)進(jìn)行奇異值分解，獲得奇異值矩陣

(10)

式中：S1、S2為非零奇異值的對角陣；U1、U2及V1、V2分別為左、右單位正交奇異陣。奇異值降序排列，1表示主要信號，2表示噪聲信號。

由于wk和vk互不相關(guān)，由離散隨機(jī)狀態(tài)模型的性質(zhì)可得

(11)

比較式(10)和式(11)可得：

(12)

再由式(11)可得

T2|a+1=OaAΓa.

(13)

將式(12)代入式(13)可得

(14)

最后根據(jù)式(2)—(4)計算低頻振蕩模態(tài)參數(shù)。

2 模態(tài)參數(shù)自動識別

2.1 基于雙協(xié)方差漢克爾矩陣的虛假模態(tài)剔除法

在不確定系統(tǒng)階數(shù)及噪聲干擾的情況下，SSI算法的模態(tài)階數(shù)n往往設(shè)置得遠(yuǎn)大于真實模態(tài)，以獲取所有的模態(tài)信息，導(dǎo)致識別結(jié)果會包含許多虛假模態(tài)。傳統(tǒng)的SSI是采用穩(wěn)定圖法，將不同階數(shù)下的分解結(jié)果繪制于同一張圖中。穩(wěn)定圖的橫坐標(biāo)為頻率，縱坐標(biāo)為模態(tài)階數(shù)，相鄰兩階次計算的模態(tài)參數(shù)差值小于設(shè)定閾值時，合并為一點，然后依靠經(jīng)驗人工選擇穩(wěn)定點，其穩(wěn)定點對應(yīng)的即為真實模態(tài)。顯然，該方法需要操作員有豐富的經(jīng)驗，操作員的能力對最終模態(tài)結(jié)果影響很大，無法實現(xiàn)真實模態(tài)自動提取。為此，本文提出了基于雙協(xié)方差的漢克爾矩陣虛假模態(tài)剔除法(false model rejection based on double covariance Hankel matrix，F(xiàn)MRDCHM)。

由前述可知，式(5)—(7)中漢克爾矩陣的行a、列b是可變量，不同的a、b取值所獲得的模態(tài)參數(shù)存在較大差異。但試驗結(jié)果表明，當(dāng)階數(shù)設(shè)置大于真實模態(tài)數(shù)時，分解結(jié)果均包含所有真實模態(tài)參數(shù)。文獻(xiàn)[16]中提到，低頻振蕩的低階模態(tài)數(shù)通常會小于5，即能量占比大的主要模態(tài)數(shù)會小于5。因此，本文將模態(tài)階數(shù)設(shè)為50，即可包含所有的低頻振蕩模態(tài)。FMRDCHM的基本原理是不同的漢克爾矩陣分解得到的結(jié)果不同，但不同階次下的虛假極點與物理極點間的距離會大于物理極點與物理極點間的距離，即

|λi,1-λi,2|<|λi,1-λj,2|.

(15)

式中：λi,1為漢克爾矩陣H1的物理極點；λi,2為漢克爾矩陣H2與之對應(yīng)的物理極點；λj,2為H2的虛假極點。通過試驗發(fā)現(xiàn)，在迭代過程中，真實模態(tài)一直存在且對其的辨識結(jié)果保持相對穩(wěn)定，而由噪聲引起的虛假模態(tài)的辨識結(jié)果波動很大；因此，以極點間的距離判定其是否為虛假極點。若極點間的距離大于1，則判定為虛假極點；反之為物理極點。

極點間的距離

(16)

式中：α為權(quán)重；W為容差；f和ξ分別代表頻率和阻尼比；下標(biāo)o、t分別代表基于H1和H2獲得的參數(shù)?？紤]到模態(tài)分析中，頻率穩(wěn)定至關(guān)重要，是系統(tǒng)穩(wěn)定的先決因素[17]，因此αf應(yīng)大于αξ，本文中αf取0.7，αξ取0.3，αf+αξ=1；另一方面，不同階次下的模態(tài)結(jié)果中，阻尼比波動遠(yuǎn)大于頻率波動，且頻率誤差較小，為反映阻尼比和頻率間的差異，Wf取值應(yīng)小于Wξ，本文中Wf、Wξ分別取0.05、0.1。為了減少計算量，可依據(jù)低頻振蕩的模態(tài)頻率和阻尼比范圍，剔除0.1～3 Hz頻率范圍外和0%～20%阻尼比范圍外的模態(tài)，再使用FMRDCHM。

FMRDCHM的流程如圖1所示，基本步驟如下：①首先根據(jù)類噪聲數(shù)據(jù)構(gòu)造2個行列不同的漢克爾矩陣H1和H2；②采用基于協(xié)方差驅(qū)動的SSI算法計算，得到2組不同階次下的模態(tài)識別結(jié)果，剔除頻率在0.1～3 Hz外的極點，剔除阻尼比在0～20%外的極點，將篩選后位于低頻振蕩模態(tài)范圍內(nèi)的2組極點，分別設(shè)定為參考組極點和驗證組極點；③根據(jù)計算同一階次下極點間距離，在驗證組中尋找與參考組對應(yīng)的最近的極點，若極點間的距離大于1，則判定為虛假極點，剔除該對極點。

2.2 基于系統(tǒng)聚類的物理模態(tài)提取

在剔除虛假模態(tài)后，進(jìn)行物理模態(tài)參數(shù)提取。本文中，為了有效確定模態(tài)階數(shù)及模態(tài)參數(shù)，采用基于系統(tǒng)聚類的物理模態(tài)提取算法。

系統(tǒng)聚類是模糊聚類中的常用方法，其基本思想是：先將每個樣本獨自看成一類，根據(jù)所定義類的距離，計算各類中的距離，將距離最小的兩類合并成一個新類，如此循環(huán)，直至合并的類滿足結(jié)束條件。本文將合并的結(jié)束條件定義為：

min|De|>1，

(17)

DA,B=min{da′,b′|a′∈A,b′∈B}.

(18)

式中：De為第e次迭代計算的距離矩陣；DA,B為類A、B間所有樣本間的最短距離；da′,b′表示類A中的樣本a′與類B中樣本b′間的距離。系統(tǒng)聚類流程如圖2所示。

基本步驟如下：

a)初始化：將FMRDCHM篩選得到的q個極點，分成q類，即C1,0,C2,0,…,Cq,0，C表示類別，下標(biāo)第1個數(shù)字表示類別數(shù)，下標(biāo)第2個數(shù)字表示當(dāng)前迭代次數(shù)；設(shè)定結(jié)束條件，令當(dāng)前迭代次數(shù)e=0。根據(jù)計算類之間的距離，得到距離矩陣

圖1 FMRDCHM流程Fig.1 Flowchart of FMRDCHM

圖2 系統(tǒng)聚類流程Fig.2 Flowchart of hierarchical clustering

(19)

求距離矩陣中最小值(對角元素除外)，假設(shè)為Ci,0與Cj,0間的距離di,j，若di,j<1，則將Ci,0與Cj,0合并為新類Cij,1，建立新類C1,1,C2,1,…,Cq-1,1，即第1次迭代得到的1，2，…，q-1類。

b)根據(jù)計算合并后新類間的距離，得到新距離矩陣D1，設(shè)D1中最小元素為DA,B，且DA,B<1，則將類A、B合并為新類。

c)令e=e+1，跳轉(zhuǎn)步驟b)，重復(fù)計算，當(dāng)min|De|>1，停止迭代。

d)統(tǒng)計最終類數(shù)及每類包含的極點數(shù)，若極點數(shù)大于閾值n/4，則判定該類為物理模態(tài)，并以該類的平均值作為最終模態(tài)參數(shù)。

3 傳遞函數(shù)生成信號

為了通過類噪聲信號評估本文方法的模態(tài)參數(shù)估計能力，首先在基于傳遞函數(shù)的模型上測試該方法，該模型產(chǎn)生的類噪聲信號包含0.5 Hz和0.8 Hz 2種模態(tài)。模態(tài)的阻尼比分別為ξ1=2%和ξ2=5%，如圖3所示。仿真模型使用均值為0、方差為1的高斯白噪聲來模擬電力系統(tǒng)中負(fù)荷隨機(jī)波動，用高斯白噪聲作為模型的激勵信號，并采用傳遞函數(shù)來模擬電力系統(tǒng)。再將2個傳遞函數(shù)的響應(yīng)信號乘以圖3中給出的系數(shù)，然后求和形成具有2種模態(tài)的類噪聲輸出信號。在該測試中，仿真時間為10 min，采樣頻率為100 Hz，生成的類噪聲信號如圖4所示，其中，s為傳遞函數(shù)的復(fù)參數(shù)，ω為模態(tài)頻率對應(yīng)的角速度，ξ為阻尼比。

圖3 基于傳遞函數(shù)的包含2種模態(tài)的測試模型Fig.3 Test model based on transfer function with two modes

圖4 基于傳遞函數(shù)生成的類噪聲信號Fig.4 Ambient signals generated by transfer function

將2個漢克爾矩陣的行數(shù)設(shè)為100和80，模型階數(shù)設(shè)為50。根據(jù)傳統(tǒng)SSI得到的原始極點圖如圖5所示。從圖5可知，虛假極點與物理極點交雜混疊，難以區(qū)分。

圖5 原始極點分布Fig.5 Pole distribution

采用傳統(tǒng)穩(wěn)定圖法進(jìn)行虛假極點剔除，以頻率為橫坐標(biāo)，模型階次為縱坐標(biāo)，建立穩(wěn)定圖，所得結(jié)果如圖6所示。由圖6可發(fā)現(xiàn)，相比于圖5的模態(tài)識別結(jié)果中夾雜許多虛假極點，傳統(tǒng)穩(wěn)定圖法所得圖像已清晰許多，大部分虛假極點被剔除，但仍有部分虛假極點存在于最終結(jié)果中，被判定為物理極點，影響最終的模態(tài)識別結(jié)果。而基于本文方法得到的穩(wěn)定圖顯然更加清晰，如圖7所示，所有虛假極點均被剔除，所有物理極點均為真實模態(tài)對應(yīng)的極點，剔除效果理想。

圖6 傳統(tǒng)方法穩(wěn)定圖Fig.6 Stabilization diagram obtained by traditional method

圖7 本文方法穩(wěn)定圖Fig.7 Stabilization diagram obtained by the proposed method

采用傳統(tǒng)穩(wěn)定圖法和本文方法進(jìn)行模態(tài)辨識，結(jié)果見表1。

表1 傳遞函數(shù)生成類噪聲信號的模態(tài)辨識結(jié)果Tab.1 Modal identification results of ambient signals generated by transfer function

基于傳統(tǒng)穩(wěn)定圖法，模態(tài)1和2的估計值都較為精確，但若篩選經(jīng)驗不足則會出現(xiàn)虛假模態(tài)，虛假模態(tài)的出現(xiàn)會對后續(xù)制訂低頻振蕩阻尼調(diào)制策略帶來不利影響。通過本文方法計算得到的估計值均接近理論值：對于頻率估計值，2種模態(tài)的頻率估計誤差均小于0.003 Hz，估計誤差幾乎可以忽略不計；而相比于頻率，阻尼比的估計誤差相對較大，但顯然也符合精度要求。更重要的是，本文方法識別結(jié)果沒有包含虛假模態(tài)，不會影響后續(xù)阻尼控制策略的制訂。即便在模態(tài)混疊的情況下，本文方法也可以精確識別低頻振蕩模態(tài)參數(shù)。

4 IEEE WSCC 3機(jī)9節(jié)點模型及算法抗噪性能測試

前一個算例是基于理想的傳遞函數(shù)模型，在此算例中，采用由IEEE WSCC 3機(jī)9節(jié)點模型生成的類噪聲信號來驗證算法的性能，IEEE WSCC 3機(jī)9節(jié)點模型會更接近真實的電力系統(tǒng)。該系統(tǒng)在電力系統(tǒng)分析工具箱(power system analysis toolbox，PSAT)中構(gòu)造，參數(shù)和結(jié)構(gòu)設(shè)置根據(jù)文獻(xiàn)[18]設(shè)置，模型如圖8所示。在電力系統(tǒng)中，類噪聲信號通常由負(fù)載中的實時隨機(jī)波動引起，為了模擬電力系統(tǒng)中負(fù)載的隨機(jī)波動，在該模型節(jié)點5、6、8處疊加30 dB的高斯白噪聲[19]。

圖8 IEEE WSCC 3機(jī)9節(jié)點系統(tǒng)Fig.8 IEEE WSCC three-machine nine-bus test system

同時基于PSAT自帶的小穩(wěn)定性分析功能進(jìn)行線性化特征值分析(linear eigenvalue analysis，LEA)，為所計算的模態(tài)參數(shù)提供可靠的參考值。但是LEA是一種基于模型的方法，當(dāng)模型較大時，計算成本高且不能滿足實時性。根據(jù)LEA結(jié)果，系統(tǒng)有2種模態(tài)：頻率f1=1.205 Hz、f2=1.831 Hz；阻尼比ξ1=2.514%、ξ2=6.125%。此外，根據(jù)PSAT產(chǎn)生的特征值報告，發(fā)電機(jī)3和發(fā)電機(jī)2各有一種模態(tài)，因此本文選擇發(fā)電機(jī)3和發(fā)電機(jī)2間的相對角速度ω32作為分析信號，如圖9所示。

圖9 IEEE WSCC 3機(jī)9節(jié)點系統(tǒng)生成的ω32信號Fig.9 Ambient signalsω32generated by IEEE WSCC three-machine nine-bus system

模態(tài)辨識結(jié)果見表2。對于頻率估計值，2種方法的真實模態(tài)1和2的結(jié)果均接近LEA結(jié)果，即頻率誤差很??；對于阻尼估計值，相比于傳統(tǒng)穩(wěn)定圖法的SSI，本文方法更接近LEA結(jié)果，即誤差更小。更重要的是，本文方法不會出現(xiàn)虛假模態(tài)3，因此本文方法的模態(tài)識別結(jié)果顯然更有優(yōu)勢。雖傳統(tǒng)穩(wěn)定圖法的平均計算時間(0.560 4 s)少于本文方法(0.616 9 s)，但犧牲較小的計算成本，獲得更為精確的結(jié)果，顯然是可以接受的。因此，本文方法更優(yōu)于傳統(tǒng)穩(wěn)定圖法SSI。

表2 不同方法下的模態(tài)辨識結(jié)果Tab.2 Modal identification resultsobtained by different methods

為了測試本文方法的抗噪能力，在測試信號上分別疊加信噪比為40 dB、35 dB、30 dB的高斯白噪聲，模態(tài)辨識結(jié)果見表3。由表3可知，雖然噪聲水平在不斷增加，但頻率的估計結(jié)果幾乎與無噪聲情況下的結(jié)果相同，阻尼比的估計結(jié)果相較于無噪聲情況下結(jié)果的最大誤差也僅為0.21%，說明本文方法具有良好的抗噪性。

5 結(jié)束語

本文針對類噪聲信號，提出了一種新的低頻振蕩模態(tài)參數(shù)估計算法?；陔p協(xié)方差的SSI算法，結(jié)合系統(tǒng)聚類方法，實現(xiàn)低頻振蕩模態(tài)辨識。本文方法無需先驗知識，可實現(xiàn)自動定階，且基于類噪聲信號，可在低阻尼或弱阻尼振蕩發(fā)生前實現(xiàn)安全預(yù)警。

表3 不同信噪比下的模態(tài)辨識結(jié)果Tab.3 Modal identification results under different SNRs

本文方法采用基于傳遞函數(shù)生成信號和基于IEEE WSCC 3機(jī)9節(jié)點系統(tǒng)的仿真信號進(jìn)行測試。仿真結(jié)果表明，該方法能從混亂無序的原始極點圖中獲取真實物理模態(tài)，相較于傳統(tǒng)SSI算法，所獲得的穩(wěn)定圖更加清晰，不包含虛假極點，且識別精度更高，識別結(jié)果僅包含真實模態(tài)，計算速度較快，具有良好的抗噪性能。