利用矩陣束算法提取水力機組參數的振蕩特性

曾保安，曾 云；張立翔（昆明理工大學建筑工程學院工程力學系，昆明 650500）

前言

水力發電機組的穩定性是衡量機組性能的指標之一，機組的穩定性將關系到整個電站及其水工廠房的安全性、國民經濟利益和整個電力系 統能否安全運行。因此，機組穩定性受到學者們廣泛關注，而表征穩定性的參數分別是振動、擺度和壓力脈動，其中振動是機組穩定運行最重要的指標。

水力發電機組的振動按其形式可分為水力振動、機械振動和電磁振動。這三種振動耦合作用于水力機組，以致造成水力機組振動機理很復雜。雖然國內外很多研究者對其做出了巨大的努力，但對其振動機理的理論研究仍然欠缺。一些研究者試圖從振動信號中提取機組的振蕩信息，找出誘發機組振蕩的誘因，以便達到改善機組穩定性的效果。目前對振動信號的提取方法有傅里葉變換法、短時傅里葉變換法、HHT變換法、小波分析法、prony算法、矩陣束算法等。其中傅里葉變換法不能反應振蕩的阻尼特性和瞬時頻率；HHT變換法會出現模態混疊的現象；小波分析方法是目前信號分析比較完備的方法，但是存在波基選取困難的局限；文獻[1]、[2]用prony方法對電力系統暫態或者大規模系統干擾進行穩定性研究，發現該方法對噪音很敏感，抗噪能力很差。文獻[3]、[4]提出一種簡化分析的方法，不依賴于對象模型可應用于從現場實測數據提取機組參數振蕩特性，可對機組任一參數的振蕩特性及其影響因素進行詳細的分析。文獻[5]、[6]研究的是矩陣束算法，此種方法最關鍵的步驟在于極點求解時進行了去噪處理，采用奇異值分解(SVD)和降低秩的方法把數據中隱含的噪音過濾掉，以免噪音產生虛假極點，其抗噪能力有很大提高。文獻[8～12]將幾種不同形式的矩陣束算法做了計算精度和計算復雜程度的比較性研究。驗證了矩陣束算法具有極強的抗噪能力和廣泛使用性。

基于以上的研究結果，本文將矩陣束算法做了一些改進，并用改進的矩陣束算法提取同一工況下水力發電機組角速度、功角、有功、進口處水頭、流量和水輪機出力等主要參數的模態信息，進而建立運動模型，初始工況、擾動強度和擾動方向是參數振蕩特性的主要影響因素，本文將對這些主要因素如何定量影響振幅、頻率和衰減因子的變化進行全面分析，總結出模態信息之間存在的耦聯關系，這些研究和探討為控制和設計提供依據。

1 矩陣束算法

1.1 矩陣束基本算法

（1）模態數的確定。通過采樣數據y(1)，y(2),…，y(N)構造一個(N-L+1)×L階 Hankel矩陣：

通過奇異值分解，Y=UΣVT，其中σi是其對角矩陣Σ對角線上的第i個奇異值。由于現場實測數據和模型存在噪聲，會產生虛假極點。可用歸一化奇異值(σi/σmax)≥β（β為閾值）來辨識出真實極點與虛假極點。把發生突變的歸一化奇異值作為初始閾值，逐漸增大或者減小閾值，選取信噪比達到最大時所對應的最大下標i作為最大模態數M，那么前M個極點為真實極點，其對應的歸一化奇異值作為閾值。

（2）構造2個Hankel矩陣。由左奇異矩陣V的前M個主奇異向量構成濾噪矩陣V′=[v1,v2,…,vM]，構造出兩個Hankel矩陣：

其中，V1、V2分別是V′刪掉第一行、最后一行而得到；Σ′的前M個奇異值與Σ相同，其余的均為0；

（3）極點與留數的求解。通過對矩陣束Y1和Y2定義，得到以下關系：

由（3）可知，當 λ≠zi時，對角陣Z0-λI的第i行不為零，即Y2-λY1的秩為M，當 λ=zi時，Z0-λI的第i行等于0；矩陣Y2-λY1的秩降為M-1，由相關理論得出，極點zi為矩陣束{Y2,Y1}的廣義特征值；即：

其中，Y1+是Y1的偽逆矩陣。

留數的求解可通過最小二乘法計算：

極點與留數求解完畢，利用(10)式求解振幅Ai、相位θi、頻率fi和衰減因子αi。

逼近函數為：

（4）擬合精度的衡量。擬合曲線與原始數據曲線之間的重合度用信噪比(SNR)來衡量；真實數據為y(i)，逼近數據為x(i)。

一般認為SNR大于20就是一種可以接受的逼近結果，該值越大，逼近效果越好。

1.2 算法改進

（1）重構數據。式(9)右邊的數據y(1)，y(2)，…，y(n)應使用消噪后的模型數據y′(1)，y′(2)，…，y′(n)，通過式子Y′=UΣ′VT求出矩陣Y′，把Y′反對角線上的數據取平均值作為模型數據。

（2）目前矩陣束算法所研究的振蕩模型通常是某個量圍繞平衡點的振蕩，其相應的相位角只分布在一、四象限，可用（10）式直接求出。但本文研究的振蕩模型屬于階躍式的振蕩，其相位角應根據留數的實部與虛部對應分布到四個象限上。

（3）影響矩陣束算法精度的因素

1）數據長度至少應取到包含全部振蕩信息。

2）采樣頻率保證大于2fmax，fmax是低頻率振蕩的最大頻率。

3）閾值是消除噪音合理與否的關鍵；一般選取突變的歸一化奇異值作為初始閾值，逐漸增大或者減小，信噪比最大時所對應的歸一化奇異值作為閾值。

4）Hankel矩陣列數L也會影響矩陣束算法的精度，一般取L/3-L/2。

2 水力機組模型

2.1 水力系統模型

水力系統采用一管多機的形式，如圖1所示，機組2和機組3假設為穩定運行狀態。鋼管與支管可用等效管i等效。

圖1 水力系統模型

其中，x=[x1x2x3x4x5]，x4=q為流量，x5=y為導葉開度。

水力系統數學模型：

fp(i)為第i路支管與鋼管的等效管道損失系數，Z(i)為第i路支管與鋼管的等效管道涌浪阻抗，T(i)為第i路支管與鋼管的等效管道彈性時間，y為主接力器位移，fpT為隧道水頭損失系數，TwT為隧道水流慣性時間系數，Ty為主接力器時間常數，qi為i路支管機組的流量，qj為除第i路支管外機組的流量。

2.2 水輪機力矩計算模型

水輪機力矩模型：

式中，mt為水輪機力矩，At為水輪機增益系數，qnl為水輪機空載流量，ht為水輪機進口處的水頭，q為水輪機流量。

2.3 發電機模型

發電機采用單機無窮大三階實用模型：

mt為輸入機械力矩；Ef為勵磁電動勢；Td′0為d軸開路暫態時間常數；Tj為機組慣性時間常數ωB=314rad/s為角速度基值；ω為機組角速度；Δω=ω－1；δ為功角；Us為無窮大系統電壓；XdΣ=Xd+XT+XL；Xd為d軸電抗；XT為變壓器電抗；XL為線路等效電抗；XqΣ=Xq+XT+XL；Xq為q軸電抗；D為阻尼系數；E′q為q軸瞬變電動勢;X′dΣ=X′d+XT+XL；X′d為d軸次暫態電抗。

3 模型計算

水力系統模型、發電機模型、水輪機出力模型、典型的并聯PID調速器和PI勵磁控制器構成完整的水力機組系統進行模型計算。

3.1 發電機功角計算模型

初始工況p=0.5p.u.,目標工況pc=1p.u.，在正常調節時，其建模的計算過程如下：

（1）按頻率10Hz進行采樣，得到功角的一組仿真數據，并利用這些數據構成一個Hankel矩陣。通過奇異值分解，選取突變的歸一化奇異值0.000048作為初始閾值，逐漸增大或減小閾值，以信噪比達到最大時所對應的下標i作為最大模態數M=7，對應的閾值β= 0 .0006。

（2）構造兩個Hankel矩陣和求解模型數據。

（3）求出極點與留數。求出的極點有3個實數和2對共軛復數，通過(10)式可知，實數極點對應的頻率為0；復數極點以共軛成對形式出現，由特征值分析理論可知，每一對復數極點對應一個振蕩模態。

（4）通過(11)式把振蕩頻率為0與不為0的運動模型分別表示出來。

功角振蕩頻率為0的運動模型：

功角振蕩頻率不為0的運動模型：

功角運動模型為δ1(t+δ2)（t）即

采用上式計算并與仿真數據進行比較，結果如圖2所示，功角運動模型曲線δ(t)與仿真曲線基本上完全重合。

圖2 p=0.5p.u., pc=1p.u.時功角的響應曲線

3.2 水輪機出力計算模型

采用與功角相同的工況和計算步驟，提取水輪機出力的運動模型。

水輪機出力振蕩頻率為0的運動模型：

水輪機出力振蕩頻率不為0的運動模型：

水輪機出力的運動模型為：

采用上式計算并與仿真數據進行比較，結果如圖3所示。水輪機出力運動模型曲線mt(t)與仿真曲線基本上完全重合。

圖3 p=0.5p.u., pc=1p.u.時水輪機出力的響應曲線

采用相同的方法研究機組主要參數的運動模型時，發現模態信息之間存在耦聯關系，即：

（1）進口處的水頭、流量和水輪機出力三個水力參數的振蕩頻率相近；功角、有功和角速度三個電氣參數的振蕩頻率相近；衰減因子α2相近；各參數的衰減因子r1相近。

（2）水力參數的擾動頻率為 0.48Hz左右；電氣參數的振動頻率為：第一基頻0.51Hz左右，第二基頻1.1Hz左右；水力參數的頻率相近，電氣參數的頻率相近，說明參數的振蕩模態只與參數類型有關。因此，在調節過程中，參數產生振蕩是機組結構本身所固有的，其頻率為固有頻率，其振蕩模態為固有振蕩模態。

通過采用相同的方法提取各參數的模態信息，并建立運動模型，且與文獻[1]建立的運動模型基本一致。因此，可將水力機組參數運動模型歸結為兩種基本形態：

1）周期衰減運動模型

其中，Ai、αi、fi、θi分別為第i個振蕩模態的振幅，衰減因子、頻率、相位角。m是周期衰減運動階數。

2）過阻尼運動模型

xz為運動終值（系統變量的平衡值），Ci、ri分別是第i個過阻尼運動模態的振幅、衰減因子，n是過阻尼運動階數。

4 模型關聯仿真

4.1 初始工況對各參數模態信息的影響

初始工況p對功角、水輪機出力模態信息的影響分別如表1、表2，其中，目標工況pc=1p.u.不變。

表1 初始工況對功角模態信息的影響

表2 初始工況對水輪機出力模態信息的影響

采用相同的方法提取其他各電氣參數、水力參數的模態信息隨初始工況變化分別與表1、表2中的變化規律一致。即：

（1）各參數的過阻尼運動的階數、周期衰減運動的階數不隨初始工況而變化，各參數的信噪比隨著初始工況減小而減小。

（2）各主要參數的振蕩頻率受初始工況影響小；電氣參數之間的頻率是相近的，水力參數之間頻率也相近。

（3）各參數的各階振幅隨初始工況減小而增大。

（4）各參數的衰減因子r1和α1隨初始工況減小而減小；各參數的衰減因子r2隨初始工況減小而增大；電氣參數的衰減因子α2隨初始工況減小而增大。

4.2 擾動強度對各參數模態信息的影響

初始工況p=0.7p.u.，擾動強度Δp對功角、水輪機出力模態信息的影響分別如表3、表4。

表3 擾動強度對功角模態信息的影響

表4 擾動強度對水輪機出力模態信息的影響

采用相同的方法提取其他各電氣參數、水力參數的模態信息隨擾動強度的變化分別與表3、表4中的變化規律一致。即：

（1）各參數的過阻尼運動的階數、周期衰減運動的階數不隨擾動強度而變化，各參數的信噪比隨著擾動強度增大而減小。

（2）各參數的振蕩頻率受擾動強度影響小，電氣參數的振蕩頻率相近，水力參數的振蕩頻率也相近，且電氣參數的振蕩頻率比水力參數的振蕩頻率多一階。

（3）各參數的振幅隨擾動強度增大而增大，不受擾動方向的影響。

（4）各參數的衰減因子r1相近；電氣參數的衰減因子α2相近；各參數的各階衰減因子隨正向擾動強度增大而減小；各階衰減因子隨負向擾動強度增大而增大。

（5）改變初始工況，Δp與表3、表4一致，可以得出各參數的模態信息隨擾動強度的變化不受初始工況影響。

5 結論

（1）改進的MP算法用于提取水力機組各參數模態信息有效性得以驗證，體現出良好的抗噪能力；該方法可通過實測數據提取各參數模態信息，不依賴于對象系統模型，因此，可用于電網和機組的在線監測和故障診斷。

（2）通過模態信息建立的兩種基本運動形態，能直觀、有效地對水力機組的振蕩特性進行描述；可通過振動頻率推斷出誘發振動的誘因。

（3）探討了各種影響因素下各參數模態信息的變化趨勢和耦聯關系，為控制設計提供依據。

[1]Hauer J F, Demeure C J, Scharf L L. Initial results in prony analysis of power system response signals [J].IEEE Trans on Power Systems, 1990, 5(1): 80-89.

[2]Grund C E, Paserba J J, J F Hauer. Comparison of prony and eigenanalysis for power system control design[J]. IEEE Transactions on Power Systems,1993,8 (3): 964-971.

[3]曾 云, 張立翔, 王煜. 水輪發電機組振蕩特性的簡化分析方法[J]. 電機與控制學報, 2009, 13(增刊1): 25-29.

[4]曾云, 沈祖詒, 曹林寧. 低頻振蕩下發電機響應特性的量化分析[J]. 電力系統及其自動化學報,2008, 20(6): 83-87.

[5]Tapan K.Sarkar , odilon Pereira. Using the matrix pencil method to estimate the parameters of a sum of complex exponentials[J]. IEEE Antennas and Propagation Magazine, 1995, 37(1): 48-55 .

[6]Tapan Kumar Sarkar, Jinhwan Koh. Application of the matrix pencil method for estimating the SEM(Singularity Expansion Method) poles of source-free transient responses from multiple look directions[J]. IEEE Transactions on Antennas and Propagation, 2000, 48(4): 612-618.

[7]Kunder P. 電力系統穩定與控制[M]. 北京:中國電力出版社, 2002.

[8]Muhammad Faisal Khan, Muhammad Tufail.Comparative analysis of various matrix pencil methods for direction of arrival estimation[C].Image Analysis and Signal Processing (IASP),International Conference on,Islamabad,Pakistan,9-11 April 2010.

[9]朱瑞可, 李興源. 矩陣束算法在同步電機參數辨識中的應用. 電力系統及其自動化學報, 2012, 36(6): 52-56.

[10]Fran?ois Sarrazin, Ala Sharaiha. Comparison between matrix pencil and prony methods applied on noisy antenna responses[C]. Loug hborough Antennas & Propagation Conference, Loughborough, UK, 14-15 November 2011.

[11]Yanhui Liu, Zaiping Nie. Reducing the number of elements in a linear antenna array by the matrix pencil method[J]. IEEE Transactions on Antennas and Propagation,2008, 56(9):2955-2962.

[12]Nuri Yilmazer, Jinhwan Koh. Utilization of a unitary transform for efficient computation in the matrix pencil method to find the direction of arrival[J]. IEEE Transactions on Antennas and Propagation, 2006, 54(1): 175-181.