基于雙協方差隨機子空間識別的類噪聲數據低頻振蕩辨識

林偉斌，季天瑤，張祿亮

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

目前，電網規模不斷擴大，系統間的互聯程度不斷增大，越來越多的高倍率快速勵磁裝置投入使用，導致系統弱阻尼低頻振蕩出現概率更高[1-2]，極大地危害電網的穩定性，也限制了互聯系統傳輸容量最大化。因此，監控和分析低頻振蕩模態參數，對確保電力系統的安全穩定具有重要意義。

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

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

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

本文方法不需人工進行定階，相比于基于單一維度下的SSI方法能更有效地區分真實模態和虛假模態。采用傳遞函數構造的合成數據以及3機9節點系統仿真數據對本文方法的有效性進行驗證，并與傳統SSI算法進行對比。

1 基于協方差的SSI算法

由環境激勵的自由度離散系統的狀態方程為[13]：

(1)

式中：A∈Rn×n為系統的狀態矩陣；C∈Rl×n為系統的輸入矩陣；xk∈Rn和yk∈Rl分別為系統k時刻的狀態量和輸出量；wk∈Rn是驅動過程中的白噪聲擾動量；vk∈Rl為測量噪聲；wk和vk均為不可測量的信號，假設它們為互不相關、均值為零的白噪聲平穩序列；n為系統階數，在本文中為低頻振蕩模態數；l為系統輸出數量。

在離散域中，對狀態矩陣A進行特征值分解：

A=ψΛψ-1.

(2)

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

(3)

(4)

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

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

首先由系統輸出量構造漢克爾矩陣：

(5)

(6)

(7)

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

根據式(5)—(7)構造托普利茨矩陣：

(8)

(9)

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

(10)

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

由于wk和vk互不相關，由離散隨機狀態模型的性質可得

(11)

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

(12)

再由式(11)可得

T2|a+1=OaAΓa.

(13)

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

(14)

最后根據式(2)—(4)計算低頻振蕩模態參數。

2 模態參數自動識別

2.1 基于雙協方差漢克爾矩陣的虛假模態剔除法

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

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

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

(15)

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

極點間的距離

(16)

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

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

2.2 基于系統聚類的物理模態提取

在剔除虛假模態后，進行物理模態參數提取。本文中，為了有效確定模態階數及模態參數，采用基于系統聚類的物理模態提取算法。

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

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′間的距離。系統聚類流程如圖2所示。

基本步驟如下：

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

圖1 FMRDCHM流程Fig.1 Flowchart of FMRDCHM

圖2 系統聚類流程Fig.2 Flowchart of hierarchical clustering

(19)

求距離矩陣中最小值(對角元素除外)，假設為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)根據計算合并后新類間的距離，得到新距離矩陣D1，設D1中最小元素為DA,B，且DA,B<1，則將類A、B合并為新類。

c)令e=e+1，跳轉步驟b)，重復計算，當min|De|>1，停止迭代。

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

3 傳遞函數生成信號

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

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

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

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

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

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

圖6 傳統方法穩定圖Fig.6 Stabilization diagram obtained by traditional method

圖7 本文方法穩定圖Fig.7 Stabilization diagram obtained by the proposed method

采用傳統穩定圖法和本文方法進行模態辨識，結果見表1。

表1 傳遞函數生成類噪聲信號的模態辨識結果Tab.1 Modal identification results of ambient signals generated by transfer function

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

4 IEEE WSCC 3機9節點模型及算法抗噪性能測試

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

圖8 IEEE WSCC 3機9節點系統Fig.8 IEEE WSCC three-machine nine-bus test system

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

圖9 IEEE WSCC 3機9節點系統生成的ω32信號Fig.9 Ambient signalsω32generated by IEEE WSCC three-machine nine-bus system

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

表2 不同方法下的模態辨識結果Tab.2 Modal identification resultsobtained by different methods

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

5 結束語

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

表3 不同信噪比下的模態辨識結果Tab.3 Modal identification results under different SNRs

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